diff --git a/DESCRIPTION b/DESCRIPTION
index 472d750..3b4e63b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: piar
 Title: Price Index Aggregation
-Version: 0.9.0.9005
+Version: 0.9.0.9006
 Authors@R: c(
     person("Steve", "Martin", role = c("aut", "cre", "cph"),
            email = "marberts@protonmail.com",
diff --git a/NEWS.md b/NEWS.md
index c2fb6f8..0fcd845 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,17 +2,21 @@
 
 ## Significant changes
 
-- Many functions from `{gpindex}` are now part of `{piar}`. This gives a somewhat
+- Many functions from `{gpindex}` are now part of this package. This gives a somewhat
 more consistent experience and allows for improvements that would be difficult with
 two separate packages. In particular:
 
     - A new function `outliers()` to find extreme price relatives.
 
-    - Functions for multilateral indexes.
+    - Functions for multilateral GEKS indexes.
 
     - Core mathematical machinery used to aggregate indexes and make product contributions.
     
     - A new vignette outlining some theory that unlies how indexes are calculated.
+    
+- The argument `r` in `elementary_index()`, `aggregate()`, `mean()`, `update()`,
+and `impute_prices()` is deprecated and will be removed in a future version.
+Use `order` instead.
 
 ## Improvements
 
diff --git a/R/aggregate.piar_index.R b/R/aggregate.piar_index.R
index 09d43d4..87ca0af 100644
--- a/R/aggregate.piar_index.R
+++ b/R/aggregate.piar_index.R
@@ -71,8 +71,8 @@
 #' @param na.rm `[logical(1)]` Should missing values be removed? By default,
 #'   missing values are not removed. Setting `na.rm = TRUE` is equivalent to
 #'   overall mean imputation.
-#' @param r `[numeric(1)]` Order of the generalized mean to aggregate index
-#'   values. 0 for a
+#' @param order,r `[numeric(1)]` Order of the generalized mean to aggregate
+#'   index values. 0 for a
 #'   geometric index (the default for making elementary indexes), 1 for an
 #'   arithmetic index (the default for aggregating elementary indexes and
 #'   averaging indexes over subperiods), or -1 for a harmonic index (usually for
@@ -149,11 +149,15 @@ aggregate.chainable_piar_index <- function(
   pias2 = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   include_ea = TRUE,
   duplicate_contrib = c("sum", "make.unique"),
   impute_rules = NULL
 ) {
+  if ("r" %in% names(sys.call())) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   chkDots(...)
   aggregate_index(
     x,
@@ -161,7 +165,7 @@ aggregate.chainable_piar_index <- function(
     pias2,
     na.rm = na.rm,
     contrib = contrib,
-    r = r,
+    order = r,
     include_ea = include_ea,
     chainable = TRUE,
     duplicate_contrib = match.arg(duplicate_contrib),
@@ -178,11 +182,15 @@ aggregate.direct_piar_index <- function(
   pias2 = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   include_ea = TRUE,
   duplicate_contrib = c("sum", "make.unique"),
   impute_rules = NULL
 ) {
+  if ("r" %in% names(sys.call())) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   chkDots(...)
   aggregate_index(
     x,
@@ -190,7 +198,7 @@ aggregate.direct_piar_index <- function(
     pias2,
     na.rm = na.rm,
     contrib = contrib,
-    r = r,
+    order = r,
     include_ea = include_ea,
     chainable = FALSE,
     duplicate_contrib = match.arg(duplicate_contrib),
@@ -206,21 +214,21 @@ aggregate_index <- function(
   pias2,
   na.rm,
   contrib,
-  r,
+  order,
   include_ea,
   chainable,
   duplicate_contrib,
   impute_rules
 ) {
   pias <- as_aggregation_structure(pias)
-  r <- as.numeric(r)
+  order <- as.numeric(order)
   has_contrib <- !is.null(x$contrib) && contrib
   res <- .aggregate(
     x,
     pias,
     na.rm,
     has_contrib,
-    r,
+    order,
     include_ea,
     chainable,
     duplicate_contrib,
@@ -246,7 +254,7 @@ aggregate_index <- function(
       pias2,
       na.rm,
       has_contrib,
-      -r,
+      -order,
       include_ea,
       chainable,
       duplicate_contrib,
@@ -259,7 +267,7 @@ aggregate_index <- function(
         res2$contrib,
         res$index,
         res2$index,
-        r = 0
+        order = 0
       )
     }
     res$index[] <- (res$index * res2$index)^0.5
@@ -279,7 +287,7 @@ aggregate_index <- function(
   pias,
   na.rm,
   has_contrib,
-  r,
+  order,
   include_ea,
   chainable,
   duplicate_contrib,
@@ -311,7 +319,7 @@ aggregate_index <- function(
       nodes <- unname(pias$child[[i - 1L]])
       rel[[i]] <- vapply(
         nodes,
-        \(z) gmean(rel[[i - 1L]][z], w[[i - 1L]][z], r, na.rm),
+        \(z) gmean(rel[[i - 1L]][z], w[[i - 1L]][z], order, na.rm),
         numeric(1L)
       )
       if (has_contrib) {
@@ -322,7 +330,7 @@ aggregate_index <- function(
               con[[i - 1L]][nodes[[j]]],
               rel[[i - 1L]][nodes[[j]]],
               w[[i - 1L]][nodes[[j]]],
-              r,
+              order,
               rel[[i]][j],
               duplicate_contrib
             )
@@ -343,7 +351,7 @@ aggregate_index <- function(
 
     # Price update weights for all periods after the first.
     if (chainable) {
-      pias$weights <- update_weights(rel[[1L]], w[[1L]], r)
+      pias$weights <- update_weights(rel[[1L]], w[[1L]], order)
     }
 
     if (!include_ea && length(rel) > 1L) {
@@ -359,8 +367,8 @@ aggregate_index <- function(
 
 #' Aggregate product contributions
 #' @noRd
-aggregate_contrib <- function(x, rel, w, r, index, duplicate_contrib) {
-  w <- transmute_weights(rel, w, r, to = 1, mean = index)
+aggregate_contrib <- function(x, rel, w, order, index, duplicate_contrib) {
+  w <- transmute_weights(rel, w, order, to = 1, mean = index)
   res <- Map(`*`, x, w)
   if (all(lengths(res) == 0L)) {
     return(numeric(0L))
@@ -385,7 +393,7 @@ aggregate_contrib <- function(x, rel, w, r, index, duplicate_contrib) {
 
 #' Aggregate product contributions for a superlative index
 #' @noRd
-super_aggregate_contrib <- function(x, y, rel1, rel2, r) {
-  w <- transmute_weights(c(rel1, rel2), order = r, to = 1)
+super_aggregate_contrib <- function(x, y, rel1, rel2, order) {
+  w <- transmute_weights(c(rel1, rel2), order = order, to = 1)
   w[1L] * x + w[2L] * y
 }
diff --git a/R/elementary_index.R b/R/elementary_index.R
index f2fca16..27fcf2d 100644
--- a/R/elementary_index.R
+++ b/R/elementary_index.R
@@ -9,7 +9,8 @@
 #' (if `contrib = TRUE`) to `x` and `weights` grouped by `ea` and `period`. That
 #' is, for every combination of elementary aggregate and time period,
 #' `elementary_index()` calculates an index based on a generalized mean of
-#' order `r` and, optionally, percent-change contributions. Product names should
+#' order `order` and, optionally, percent-change contributions. Product names
+#' should
 #' be unique within each elementary aggregate at each time period when making
 #' contributions and, if not, are
 #' passed to [make.unique()] with a warning. The default
@@ -81,8 +82,8 @@
 #'   missing values
 #'   are not removed. Setting `na.rm = TRUE` is equivalent to overall mean
 #'   imputation.
-#' @param r `[numeric(1)]` Order of the generalized mean to aggregate price
-#'   relatives. 0 for a
+#' @param order,r `[numeric(1)]` Order of the generalized mean to aggregate
+#'   price relatives. 0 for a
 #'   geometric index (the default for making elementary indexes), 1 for an
 #'   arithmetic index (the default for aggregating elementary indexes and
 #'   averaging indexes over subperiods), or -1 for a harmonic index (usually for
@@ -163,6 +164,9 @@
 #'   r = 1
 #' )
 elementary_index <- function(x, ...) {
+  if ("r" %in% ...names()) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   UseMethod("elementary_index")
 }
 
@@ -184,7 +188,8 @@ elementary_index.numeric <- function(
   chainable = TRUE,
   na.rm = FALSE,
   contrib = FALSE,
-  r = 0
+  order = 0,
+  r = order
 ) {
   chkDots(...)
   if (!is.null(weights)) {
diff --git a/R/emean.R b/R/emean.R
index 21f66bd..4923356 100644
--- a/R/emean.R
+++ b/R/emean.R
@@ -2,6 +2,9 @@
 #'
 #' Calculate the component-wise extended mean.
 #'
+#' The extended mean is also called the difference mean, Stolarsky mean, or
+#' extended mean-value mean; see Bullen (2003, p. 393) for details.
+#'
 #' Both `x` and `y` should be strictly positive. This is not enforced, but the
 #' results may not make sense when the extended mean is not defined. The usual
 #' recycling rules apply when `x` and `y` are not the same length.
@@ -74,19 +77,20 @@ emean <- function(
   # Recycling x and y here avoids multiple warnings if one is not a multiple
   # length of the other.
   if (length(x) > length(y)) {
-    if (length(x) %% length(y) != 0) {
-      warning("length of `x` is not a multiple of length of `y")
-      y <- rep_len(length(x))
+    if (length(y) > 0 && length(x) %% length(y) != 0) {
+      warning("length of `y` is not a multiple of length of `x`")
+      y <- rep_len(y, length(x))
     }
   } else if (length(x) < length(y)) {
-    if (length(y) %% length(x) != 0) {
-      warning("length of `y` is not a multiple of length of `x")
-      x <- rep_len(length(y))
+    if (length(x) > 0 && length(y) %% length(x) != 0) {
+      warning("length of `x` is not a multiple of length of `y`")
+      x <- rep_len(x, length(y))
     }
   }
 
-  if (length(tol) > max(length(x), length(y))) {
-    stop("'tol' cannot be longer than 'x' or 'y'")
+  max_len <- max(length(x), length(y))
+  if (max_len > 0 && length(tol) > max_len) {
+    stop("`tol` cannot be longer than `x` or `y`")
   }
 
   if (r == 0 && s == 0) {
diff --git a/R/gmean.R b/R/gmean.R
index 8c9ab32..1a67f28 100644
--- a/R/gmean.R
+++ b/R/gmean.R
@@ -31,8 +31,8 @@
 #' @param x `[numeric > 0]` A strictly positive numeric vector.
 #' @param weights `[numeric >= 0]` A positive numeric vector of weights, the
 #'   same length as `x`. The default is to equally weight each element of `x`.
-#' @param order `[numeric(1)]` A finite number giving the order of the
-#'   generalized mean. The default calculates an arithmetic mean.
+#' @param order `[numeric(1)]` A finite number giving the order (or exponent) of
+#'   the generalized mean. The default calculates an arithmetic mean.
 #' @param na.rm `[logical(1)]` Should missing values be removed? By default,
 #'   missing values are not removed.
 #'
diff --git a/R/impute_prices.R b/R/impute_prices.R
index ec22030..0de56e0 100644
--- a/R/impute_prices.R
+++ b/R/impute_prices.R
@@ -59,8 +59,8 @@
 #'   product weights), or something that can be coerced into one. The default is
 #'   to give each price equal weight. This is evaluated in `x` for the data
 #'   frame method.
-#' @param r `[numeric(2)]` A pair of numeric values. The first gives the order
-#'   of the generalized-mean price index used to calculate the
+#' @param order,r `[numeric(2)]` A pair of numeric values. The first gives the
+#'   order of the generalized-mean price index used to calculate the
 #'   elementary price indexes, defaulting to a geometric index. The second
 #'   gives the order of the generalized-mean price index used to aggregate the
 #'   elementary price indexes, defaulting to an arithmetic index. Other values
@@ -123,6 +123,9 @@
 #'   method = "overall-mean"
 #' )
 impute_prices <- function(x, ...) {
+  if ("r" %in% ...names()) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   UseMethod("impute_prices")
 }
 
@@ -142,7 +145,8 @@ impute_prices.matrix <- function(
   ea = NULL,
   weights = NULL,
   pias = NULL,
-  r = c(0, 1),
+  order = c(0, 1),
+  r = order,
   method = c("overall-mean", "carry-forward"),
   impute_rules = NULL
 ) {
@@ -192,7 +196,7 @@ impute_prices.matrix <- function(
         ea = ea[[t]],
         weights = weights[[t]],
         na.rm = TRUE,
-        r = r[1L]
+        order = r[1L]
       )
       time(index) <- names(res)[t]
       if (!is.null(pias)) {
@@ -200,10 +204,10 @@ impute_prices.matrix <- function(
           index,
           pias,
           na.rm = TRUE,
-          r = r[2L],
+          order = r[2L],
           impute_rules = impute_rules
         )
-        pias <- update(pias, index, r = r[2L])
+        pias <- update(pias, index, order = r[2L])
       }
       eas <- if (!is.null(ea)) {
         match(as.character(ea[[t]][impute]), index$levels)
@@ -240,7 +244,8 @@ impute_prices.numeric <- function(
   ea = NULL,
   weights = NULL,
   pias = NULL,
-  r = c(0, 1),
+  order = c(0, 1),
+  r = order,
   method = c("overall-mean", "carry-forward", "carry-backward"),
   impute_rules = NULL
 ) {
@@ -296,7 +301,7 @@ impute_prices.numeric <- function(
         ea = ea[[t]],
         weights = weights[[t]],
         na.rm = TRUE,
-        r = r[1L]
+        order = r[1L]
       )
       time(index) <- names(res)[t]
       if (!is.null(pias)) {
@@ -304,10 +309,10 @@ impute_prices.numeric <- function(
           index,
           pias,
           na.rm = TRUE,
-          r = r[2L],
+          order = r[2L],
           impute_rules = impute_rules
         )
-        pias <- update(pias, index, r = r[2L])
+        pias <- update(pias, index, order = r[2L])
       }
       eas <- if (!is.null(ea)) {
         match(as.character(ea[[t]][impute]), index$levels)
diff --git a/R/mean.piar_index.R b/R/mean.piar_index.R
index f92c4a3..190aa2a 100644
--- a/R/mean.piar_index.R
+++ b/R/mean.piar_index.R
@@ -39,8 +39,8 @@
 #'   missing values
 #'   are not removed. Setting `na.rm = TRUE` is equivalent to overall mean
 #'   imputation.
-#' @param r `[numeric(1)]` Order of the generalized mean to aggregate index
-#'   values. 0 for a
+#' @param order,r `[numeric(1)]` Order of the generalized mean to aggregate
+#'   index values. 0 for a
 #'   geometric index (the default for making elementary indexes), 1 for an
 #'   arithmetic index (the default for aggregating elementary indexes and
 #'   averaging indexes over subperiods), or -1 for a harmonic index (usually for
@@ -74,17 +74,21 @@ mean.chainable_piar_index <- function(
   window = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   duplicate_contrib = c("sum", "make.unique")
 ) {
   chkDots(...)
+  if ("r" %in% names(sys.call())) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   mean_index(
     x,
     weights,
     window = window,
     na.rm = na.rm,
     contrib = contrib,
-    r = r,
+    order = r,
     chainable = TRUE,
     duplicate_contrib = match.arg(duplicate_contrib)
   )
@@ -99,17 +103,21 @@ mean.direct_piar_index <- function(
   window = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   duplicate_contrib = c("sum", "make.unique")
 ) {
   chkDots(...)
+  if ("r" %in% names(sys.call())) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   mean_index(
     x,
     weights,
     window = window,
     na.rm = na.rm,
     contrib = contrib,
-    r = r,
+    order = r,
     chainable = FALSE,
     duplicate_contrib = match.arg(duplicate_contrib)
   )
@@ -123,7 +131,7 @@ mean_index <- function(
   window,
   na.rm,
   contrib,
-  r,
+  order,
   chainable,
   duplicate_contrib
 ) {
@@ -138,7 +146,7 @@ mean_index <- function(
     dim(weights) <- c(nlevels(x), ntime(x))
   }
 
-  r <- as.numeric(r)
+  order <- as.numeric(order)
   window <- as.integer(window %||% ntime(x))
   if (window < 1L) {
     stop("`window` must be a positive integer")
@@ -166,7 +174,7 @@ mean_index <- function(
     if (!is.null(weights)) {
       w <- split_rows(weights[, j, drop = FALSE], rows)
     }
-    res[[i]] <- mapply(gmean, rel, w, r, na.rm = na.rm, USE.NAMES = FALSE)
+    res[[i]] <- mapply(gmean, rel, w, order, na.rm = na.rm, USE.NAMES = FALSE)
     if (has_contrib) {
       con <- split_rows(x$contrib[, j, drop = FALSE], rows)
       contrib[[i]] <- mapply(
@@ -174,7 +182,7 @@ mean_index <- function(
         con,
         rel,
         w,
-        r,
+        order,
         res[[i]],
         duplicate_contrib,
         SIMPLIFY = FALSE,
diff --git a/R/update.aggregation_structure.R b/R/update.aggregation_structure.R
index 41854fd..43ba1a7 100644
--- a/R/update.aggregation_structure.R
+++ b/R/update.aggregation_structure.R
@@ -13,8 +13,8 @@
 #'   [`aggregate()`][aggregate.piar_index].
 #' @param period `[character(1)]` The time period used to price update the
 #'   weights. The default uses the last period in `index`.
-#' @param r `[numeric(1)]` Order of the generalized mean to update the weights.
-#'   The default is 1 for an arithmetic index.
+#' @param order,r `[numeric(1)]` Order of the generalized mean to update the
+#'   weights. The default is 1 for an arithmetic index.
 #' @param ... Not currently used.
 #'
 #' @returns
@@ -53,8 +53,12 @@ update.piar_aggregation_structure <- function(
   index,
   ...,
   period = NULL,
-  r = 1
+  order = 1,
+  r = order
 ) {
+  if ("r" %in% names(sys.call())) {
+    warning("`r` is deprecated and will be removed; use `order` instead")
+  }
   chkDots(...)
   index <- as_index(index, chainable = FALSE)
   r <- as.numeric(r)
@@ -70,7 +74,7 @@ update.piar_aggregation_structure <- function(
   weights(object) <- update_weights(
     index$index[, period][eas],
     object$weights,
-    r
+    order = r
   )
   object
 }
diff --git a/inst/tinytest/test-aggregate.R b/inst/tinytest/test-aggregate.R
index 372a116..1f35635 100644
--- a/inst/tinytest/test-aggregate.R
+++ b/inst/tinytest/test-aggregate.R
@@ -224,7 +224,7 @@ local({
     price_relative(price, period = period, product = product) ~
       period + business,
     contrib = TRUE,
-    r = 0.2
+    order = 0.2
   )
 
   ms_pias <- with(
@@ -235,7 +235,7 @@ local({
     )
   )
 
-  ms_index <- aggregate(ms_epr, ms_pias, r = -1.7, na.rm = TRUE)
+  ms_index <- aggregate(ms_epr, ms_pias, order = -1.7, na.rm = TRUE)
 
   res <- c(
     1,
@@ -283,11 +283,14 @@ local({
   expect_equal(as.matrix(ms_index), res)
 
   expect_equal(
-    as.matrix(aggregate(ms_index, ms_pias, r = -1.7, na.rm = TRUE)),
+    as.matrix(aggregate(ms_index, ms_pias, order = -1.7, na.rm = TRUE)),
     as.matrix(ms_index)
   )
 
-  expect_equal(aggregate(chain(ms_index), ms_pias, r = -1.7), chain(ms_index))
+  expect_equal(
+    aggregate(chain(ms_index), ms_pias, order = -1.7),
+    chain(ms_index)
+  )
 
   expect_equal(
     as.matrix(ms_index)[1, ],
@@ -295,10 +298,11 @@ local({
   )
   expect_equal(
     as.matrix(ms_index)[1, ],
-    colSums(contrib(aggregate(ms_index, ms_pias, r = -1.7)), na.rm = TRUE) + 1
+    colSums(contrib(aggregate(ms_index, ms_pias, order = -1.7)), na.rm = TRUE) +
+      1
   )
   expect_equal(
-    aggregate(ms_index, ms_pias, r = -1.7, duplicate_contrib = "sum"),
+    aggregate(ms_index, ms_pias, order = -1.7, duplicate_contrib = "sum"),
     ms_index
   )
 
@@ -320,9 +324,9 @@ local({
     as.matrix(chain(ms_index[1:3, ]))
   )
 
-  ms_index <- aggregate(ms_epr, ms_pias, r = -1.7)
+  ms_index <- aggregate(ms_epr, ms_pias, order = -1.7)
 
-  expect_equal(aggregate(ms_index, ms_pias, r = -1.7), ms_index)
+  expect_equal(aggregate(ms_index, ms_pias, order = -1.7), ms_index)
 
   expect_equal(as.matrix(ms_index)[1, ], colSums(contrib(ms_index)) + 1)
   expect_equal(
@@ -476,8 +480,8 @@ local({
 
   # Should work for a non-arithmetic index
   expect_equal(
-    chain(aggregate(epr_pop, pias, r = 3)),
-    aggregate(epr_fx, pias, r = 3)
+    chain(aggregate(epr_pop, pias, order = 3)),
+    aggregate(epr_fx, pias, order = 3)
   )
 
   # Consistency in aggregation holds with a change in base period
@@ -519,8 +523,8 @@ local({
     1:3
   )
 
-  arithmetic_contributions <- function(x, w, r = 1) {
-    (x - 1) * transmute_weights(x, w, r, to = 1)
+  arithmetic_contributions <- function(x, w, order = 1) {
+    (x - 1) * transmute_weights(x, w, order, to = 1)
   }
   index <- aggregate(epr, pias)
   expect_equal(
@@ -548,7 +552,7 @@ local({
     ea = gl(2, 2),
     period = gl(1, 4),
     contrib = TRUE,
-    r = 1
+    order = 1
   )
   epr2 <- epr1
   levels(epr2) <- 3:4
@@ -621,9 +625,9 @@ local({
   )
   levels(epr2) <- c(211, 221, 222)
 
-  index <- aggregate(merge(epr1, epr2), pias, r = 0.5)
-  index1 <- aggregate(epr1, pias, r = 0.5)
-  index2 <- aggregate(epr2, pias, r = 0.5)
+  index <- aggregate(merge(epr1, epr2), pias, order = 0.5)
+  index1 <- aggregate(epr1, pias, order = 0.5)
+  index2 <- aggregate(epr2, pias, order = 0.5)
 
   expect_equal(index[1], index1[1])
   expect_equal(index[2], index2[2])
@@ -691,11 +695,14 @@ local({
     )
   )
 
-  ms_index <- chain(aggregate(ms_epr, ms_pias, na.rm = TRUE, r = 2))
+  ms_index <- chain(aggregate(ms_epr, ms_pias, na.rm = TRUE, order = 2))
 
   ms_index2 <- unchain(ms_index[, -3])
 
-  expect_equal(chain(aggregate(ms_index2, ms_pias, r = 2))[, 3], ms_index[, 4])
+  expect_equal(
+    chain(aggregate(ms_index2, ms_pias, order = 2))[, 3],
+    ms_index[, 4]
+  )
 })
 
 # Skipping eas works.
@@ -778,7 +785,7 @@ local({
   )
   res <- aggregate(epr, pias, pias2 = pias2, na.rm = TRUE)
   res1 <- aggregate(epr, pias, na.rm = TRUE)
-  res2 <- aggregate(epr, pias2, r = -1, na.rm = TRUE)
+  res2 <- aggregate(epr, pias2, order = -1, na.rm = TRUE)
 
   expect_equal(
     sqrt(as.matrix(res1) * as.matrix(res2)),
@@ -812,8 +819,8 @@ local({
 
   # Tornqvist.
   expect_equal(
-    aggregate(epr, pias, r = 0, na.rm = TRUE),
-    aggregate(epr, pias, pias2 = pias, r = 0, na.rm = TRUE)
+    aggregate(epr, pias, order = 0, na.rm = TRUE),
+    aggregate(epr, pias, pias2 = pias, order = 0, na.rm = TRUE)
   )
 })
 
diff --git a/inst/tinytest/test-elemental_index.R b/inst/tinytest/test-elemental_index.R
index f8a4628..d51d477 100644
--- a/inst/tinytest/test-elemental_index.R
+++ b/inst/tinytest/test-elemental_index.R
@@ -7,7 +7,7 @@ epr1 <- elementary_index(ms_prices, rel ~ period + business, contrib = TRUE)
 epr2 <- elementary_index(
   ms_prices,
   rel ~ period + business,
-  r = -1,
+  order = -1,
   contrib = TRUE,
   na.rm = TRUE
 )
@@ -58,13 +58,13 @@ local({
     ms_prices2,
     rel ~ period + business,
     weights = w1,
-    r = 1
+    order = 1
   )
   p <- elementary_index(
     ms_prices2,
     rel ~ period + business,
     weights = w2,
-    r = -1
+    order = -1
   )
   expect_equal(sqrt(as.matrix(l) * as.matrix(p)), as.matrix(epr3))
 
@@ -80,12 +80,12 @@ local({
 
   sepr <- elementary_index(ms_prices2, rel ~ period + business, weights = w)
 
-  l <- elementary_index(ms_prices2, rel ~ period + business, r = 1.5)
+  l <- elementary_index(ms_prices2, rel ~ period + business, order = 1.5)
   p <- elementary_index(
     ms_prices2,
     rel ~ period + business,
     weights = w2,
-    r = -1.5
+    order = -1.5
   )
   expect_equal(sqrt(as.matrix(l) * as.matrix(p)), as.matrix(sepr))
 })
@@ -117,7 +117,7 @@ local({
   expect_error(elementary_index(1:3, period = factor(1:3, levels = numeric(0))))
   expect_error(elementary_index(1:3, ea = factor(1:3, levels = numeric(0))))
   expect_error(elementary_index(setNames(1:3, c("", 1, 2)), contrib = TRUE))
-  expect_error(elementary_index(-1:1, period = 1:3, ea = 1:3, r = 1))
+  expect_error(elementary_index(-1:1, period = 1:3, ea = 1:3, order = 1))
   expect_warning(elementary_index(
     setNames(1:3, rep(1, 3)),
     period = gl(1, 3),
diff --git a/inst/tinytest/test-emean.R b/inst/tinytest/test-emean.R
index afd8ff1..612b0ea 100644
--- a/inst/tinytest/test-emean.R
+++ b/inst/tinytest/test-emean.R
@@ -1,4 +1,4 @@
-# Extended mean satisfied key identities.
+# Extended mean satisfies key identities.
 a <- 1:5
 b <- c(1, 5, 4, 1, 2)
 
@@ -20,6 +20,7 @@ expect_equal(
   emean(a, b, order = c(1, 1)),
   emean(b, a, order = c(1, 1))
 )
+
 # Identities.
 expect_equal(
   emean(a, b, order = c(-1, 1)),
@@ -66,4 +67,30 @@ expect_equal(
 )
 
 # Errors when expected.
-expect_error(emean(a, b, tol = 1:10))
+expect_error(emean(a, b, tol = 1:10), "`tol` cannot be longer than `x` or `y`")
+expect_error(
+  emean(a, b, order = 1:3),
+  "`order` must be a pair of finite numbers"
+)
+
+# Recycling works
+expect_warning(
+  emean(1:3, 1:5),
+  "length of `x` is not a multiple of length of `y`"
+)
+expect_equal(
+  suppressWarnings(emean(1:3, 1:5)),
+  emean(c(1:3, 1:2), 1:5)
+)
+expect_equal(emean(1:5, numeric(0)), numeric(0))
+
+expect_warning(
+  emean(1:5, 1:3),
+  "length of `y` is not a multiple of length of `x`"
+)
+expect_equal(
+  suppressWarnings(emean(1:5, 1:3)),
+  emean(1:5, c(1:3, 1:2))
+)
+expect_equal(emean(numeric(0), 1:5), numeric(0))
+expect_equal(emean(numeric(0), numeric(0)), numeric(0))
diff --git a/inst/tinytest/test-mean.R b/inst/tinytest/test-mean.R
index 0285c80..ee89062 100644
--- a/inst/tinytest/test-mean.R
+++ b/inst/tinytest/test-mean.R
@@ -63,7 +63,7 @@ local({
     ms_epr,
     weights = w,
     window = 2,
-    r = 2.5,
+    order = 2.5,
     na.rm = TRUE,
     duplicate_contrib = "sum"
   )
@@ -78,7 +78,7 @@ local({
 local({
   x <- as_index(matrix(1:9, 3))
   expect_equal(mean(x, window = 1), x)
-  expect_equal(mean(x, weights = 9:1, window = 1, r = 0), x)
+  expect_equal(mean(x, weights = 9:1, window = 1, order = 0), x)
   expect_equal(mean(window(x, end(x))), window(x, end(x)))
 })
 
diff --git a/inst/tinytest/test-stack.R b/inst/tinytest/test-stack.R
index 04a4bbf..b88a3e4 100644
--- a/inst/tinytest/test-stack.R
+++ b/inst/tinytest/test-stack.R
@@ -23,7 +23,7 @@ local({
   expect_error(suppressWarnings(stack(index1, epr1)))
   time(epr2) <- time(epr1)
   expect_error(stack(epr1, epr2))
-  expect_error(stack(index1, aggregate(epr2, pias2, r = 0)))
+  expect_error(stack(index1, aggregate(epr2, pias2, order = 0)))
 })
 
 # Stacking and unstacking are opposite operations.
diff --git a/inst/tinytest/test-update.R b/inst/tinytest/test-update.R
index 4bee4a9..ff756d4 100644
--- a/inst/tinytest/test-update.R
+++ b/inst/tinytest/test-update.R
@@ -32,7 +32,7 @@ local({
     )
   )
   expect_equal(
-    weights(update(agg, index, r = 0), ea_only = FALSE),
+    weights(update(agg, index, order = 0), ea_only = FALSE),
     list(
       level1 = c("1" = NA_real_, "2" = NA),
       level2 = c("11" = NA_real_, "21" = NA, "12" = 1),
diff --git a/man/aggregate.piar_index.Rd b/man/aggregate.piar_index.Rd
index e0ee447..f822023 100644
--- a/man/aggregate.piar_index.Rd
+++ b/man/aggregate.piar_index.Rd
@@ -13,7 +13,8 @@
   pias2 = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   include_ea = TRUE,
   duplicate_contrib = c("sum", "make.unique"),
   impute_rules = NULL
@@ -26,7 +27,8 @@
   pias2 = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   include_ea = TRUE,
   duplicate_contrib = c("sum", "make.unique"),
   impute_rules = NULL
@@ -52,8 +54,8 @@ overall mean imputation.}
 \item{contrib}{\verb{[logical(1)]} Aggregate percent-change contributions in \code{x}?
 By default contributions are aggregated.}
 
-\item{r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate index
-values. 0 for a
+\item{order, r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate
+index values. 0 for a
 geometric index (the default for making elementary indexes), 1 for an
 arithmetic index (the default for aggregating elementary indexes and
 averaging indexes over subperiods), or -1 for a harmonic index (usually for
diff --git a/man/elementary_index.Rd b/man/elementary_index.Rd
index 1cac3f9..84e335f 100644
--- a/man/elementary_index.Rd
+++ b/man/elementary_index.Rd
@@ -22,7 +22,8 @@ elementary_index(x, ...)
   chainable = TRUE,
   na.rm = FALSE,
   contrib = FALSE,
-  r = 0
+  order = 0,
+  r = order
 )
 
 \method{elementary_index}{data.frame}(x, formula, ..., weights = NULL, product = NULL)
@@ -72,8 +73,8 @@ imputation.}
 \item{contrib}{\verb{[logical(1)]} Should percent-change contributions be
 calculated? The default does not calculate contributions.}
 
-\item{r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate price
-relatives. 0 for a
+\item{order, r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate
+price relatives. 0 for a
 geometric index (the default for making elementary indexes), 1 for an
 arithmetic index (the default for aggregating elementary indexes and
 averaging indexes over subperiods), or -1 for a harmonic index (usually for
@@ -102,7 +103,8 @@ wrapper that applies \code{\link[=gmean]{gmean()}} and \code{\link[=transmute_we
 (if \code{contrib = TRUE}) to \code{x} and \code{weights} grouped by \code{ea} and \code{period}. That
 is, for every combination of elementary aggregate and time period,
 \code{elementary_index()} calculates an index based on a generalized mean of
-order \code{r} and, optionally, percent-change contributions. Product names should
+order \code{order} and, optionally, percent-change contributions. Product names
+should
 be unique within each elementary aggregate at each time period when making
 contributions and, if not, are
 passed to \code{\link[=make.unique]{make.unique()}} with a warning. The default
diff --git a/man/emean.Rd b/man/emean.Rd
index 9b21dc7..196c890 100644
--- a/man/emean.Rd
+++ b/man/emean.Rd
@@ -26,6 +26,9 @@ of \code{x} and \code{y}.
 Calculate the component-wise extended mean.
 }
 \details{
+The extended mean is also called the difference mean, Stolarsky mean, or
+extended mean-value mean; see Bullen (2003, p. 393) for details.
+
 Both \code{x} and \code{y} should be strictly positive. This is not enforced, but the
 results may not make sense when the extended mean is not defined. The usual
 recycling rules apply when \code{x} and \code{y} are not the same length.
diff --git a/man/gmean.Rd b/man/gmean.Rd
index 6082f1c..434c3a1 100644
--- a/man/gmean.Rd
+++ b/man/gmean.Rd
@@ -12,8 +12,8 @@ gmean(x, weights = NULL, order = 1, na.rm = FALSE)
 \item{weights}{\verb{[numeric >= 0]} A positive numeric vector of weights, the
 same length as \code{x}. The default is to equally weight each element of \code{x}.}
 
-\item{order}{\verb{[numeric(1)]} A finite number giving the order of the
-generalized mean. The default calculates an arithmetic mean.}
+\item{order}{\verb{[numeric(1)]} A finite number giving the order (or exponent) of
+the generalized mean. The default calculates an arithmetic mean.}
 
 \item{na.rm}{\verb{[logical(1)]} Should missing values be removed? By default,
 missing values are not removed.}
diff --git a/man/impute_prices.Rd b/man/impute_prices.Rd
index b44c020..f6b82e5 100644
--- a/man/impute_prices.Rd
+++ b/man/impute_prices.Rd
@@ -23,7 +23,8 @@ impute_prices(x, ...)
   ea = NULL,
   weights = NULL,
   pias = NULL,
-  r = c(0, 1),
+  order = c(0, 1),
+  r = order,
   method = c("overall-mean", "carry-forward"),
   impute_rules = NULL
 )
@@ -36,7 +37,8 @@ impute_prices(x, ...)
   ea = NULL,
   weights = NULL,
   pias = NULL,
-  r = c(0, 1),
+  order = c(0, 1),
+  r = order,
   method = c("overall-mean", "carry-forward", "carry-backward"),
   impute_rules = NULL
 )
@@ -83,8 +85,8 @@ structure, or something that can be
 coerced into one, as made with \code{\link[=aggregation_structure]{aggregation_structure()}}. The default
 imputes from elementary indexes only (i.e., not recursively).}
 
-\item{r}{\verb{[numeric(2)]} A pair of numeric values. The first gives the order
-of the generalized-mean price index used to calculate the
+\item{order, r}{\verb{[numeric(2)]} A pair of numeric values. The first gives the
+order of the generalized-mean price index used to calculate the
 elementary price indexes, defaulting to a geometric index. The second
 gives the order of the generalized-mean price index used to aggregate the
 elementary price indexes, defaulting to an arithmetic index. Other values
diff --git a/man/mean.piar_index.Rd b/man/mean.piar_index.Rd
index d549f4e..2d37c0e 100644
--- a/man/mean.piar_index.Rd
+++ b/man/mean.piar_index.Rd
@@ -13,7 +13,8 @@
   window = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   duplicate_contrib = c("sum", "make.unique")
 )
 
@@ -24,7 +25,8 @@
   window = NULL,
   na.rm = FALSE,
   contrib = TRUE,
-  r = 1,
+  order = 1,
+  r = order,
   duplicate_contrib = c("sum", "make.unique")
 )
 }
@@ -54,8 +56,8 @@ imputation.}
 \item{contrib}{\verb{[logical(1)]} Aggregate percent-change contributions in \code{x}?
 By default contributions are aggregated.}
 
-\item{r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate index
-values. 0 for a
+\item{order, r}{\verb{[numeric(1)]} Order of the generalized mean to aggregate
+index values. 0 for a
 geometric index (the default for making elementary indexes), 1 for an
 arithmetic index (the default for aggregating elementary indexes and
 averaging indexes over subperiods), or -1 for a harmonic index (usually for
diff --git a/man/update.piar_aggregation_structure.Rd b/man/update.piar_aggregation_structure.Rd
index beb37dd..013cdad 100644
--- a/man/update.piar_aggregation_structure.Rd
+++ b/man/update.piar_aggregation_structure.Rd
@@ -4,7 +4,7 @@
 \alias{update.piar_aggregation_structure}
 \title{Update an aggregation structure}
 \usage{
-\method{update}{piar_aggregation_structure}(object, index, ..., period = NULL, r = 1)
+\method{update}{piar_aggregation_structure}(object, index, ..., period = NULL, order = 1, r = order)
 }
 \arguments{
 \item{object}{\verb{[piar_aggregation_structure]} A price index aggregation
@@ -19,8 +19,8 @@ that can be coerced into one. Usually an aggregate price index as made by
 \item{period}{\verb{[character(1)]} The time period used to price update the
 weights. The default uses the last period in \code{index}.}
 
-\item{r}{\verb{[numeric(1)]} Order of the generalized mean to update the weights.
-The default is 1 for an arithmetic index.}
+\item{order, r}{\verb{[numeric(1)]} Order of the generalized mean to update the
+weights. The default is 1 for an arithmetic index.}
 }
 \value{
 A copy of \code{object} with price-updated weights using the index
diff --git a/man/update_weights.Rd b/man/update_weights.Rd
index dfb0145..013a512 100644
--- a/man/update_weights.Rd
+++ b/man/update_weights.Rd
@@ -12,8 +12,8 @@ update_weights(x, weights = NULL, order = 1)
 \item{weights}{\verb{[numeric >= 0]} A positive numeric vector of weights, the
 same length as \code{x}. The default is to equally weight each element of \code{x}.}
 
-\item{order}{\verb{[numeric(1)]} A finite number giving the order of the
-generalized mean. The default calculates an arithmetic mean.}
+\item{order}{\verb{[numeric(1)]} A finite number giving the order (or exponent) of
+the generalized mean. The default calculates an arithmetic mean.}
 }
 \value{
 A numeric vector the same length as \code{x}.
diff --git a/tests/test-making-price-indexes.R b/tests/test-making-price-indexes.R
index c60397b..e47a286 100644
--- a/tests/test-making-price-indexes.R
+++ b/tests/test-making-price-indexes.R
@@ -132,10 +132,20 @@ ms_elementary2
 aggregate(ms_elementary2, pias, na.rm = TRUE)
 
 ## -----------------------------------------------------------------------------
-elementary_index(ms_prices, relatives ~ period + business, na.rm = TRUE, r = 1)
+elementary_index(
+  ms_prices,
+  relatives ~ period + business,
+  na.rm = TRUE,
+  order = 1
+)
 
 ## -----------------------------------------------------------------------------
-elementary_index(ms_prices, relatives ~ period + business, na.rm = TRUE, r = -1)
+elementary_index(
+  ms_prices,
+  relatives ~ period + business,
+  na.rm = TRUE,
+  order = -1
+)
 
 ## -----------------------------------------------------------------------------
 ms_prices2 <- transform(ms_prices, quantity = 10 - price)
@@ -294,7 +304,7 @@ pias <- with(
 ## -----------------------------------------------------------------------------
 paasche <- Reduce(
   stack,
-  Map(aggregate, ms_elementary, pias, na.rm = TRUE, r = -1)
+  Map(aggregate, ms_elementary, pias, na.rm = TRUE, order = -1)
 )
 
 paasche
diff --git a/tests/test-making-price-indexes.Rout.save b/tests/test-making-price-indexes.Rout.save
index 68e91a5..b487c1f 100644
--- a/tests/test-making-price-indexes.Rout.save
+++ b/tests/test-making-price-indexes.Rout.save
@@ -299,7 +299,12 @@ levels 202001    202002    202003   202004
     B5      1 1.1721550 1.0400686 4.576286
 > 
 > ## -----------------------------------------------------------------------------
-> elementary_index(ms_prices, relatives ~ period + business, na.rm = TRUE, r = 1)
+> elementary_index(
++   ms_prices,
++   relatives ~ period + business,
++   na.rm = TRUE,
++   order = 1
++ )
 Period-over-period price index for 4 levels over 4 time periods 
       time
 levels 202001     202002    202003   202004
@@ -309,7 +314,12 @@ levels 202001     202002    202003   202004
     B4    NaN        NaN       NaN 9.368610
 > 
 > ## -----------------------------------------------------------------------------
-> elementary_index(ms_prices, relatives ~ period + business, na.rm = TRUE, r = -1)
+> elementary_index(
++  ms_prices,
++  relatives ~ period + business,
++  na.rm = TRUE,
++  order = -1
++ )
 Period-over-period price index for 4 levels over 4 time periods 
       time
 levels 202001    202002    202003    202004
@@ -636,7 +646,7 @@ levels   202004
 > ## -----------------------------------------------------------------------------
 > paasche <- Reduce(
 +   stack,
-+   Map(aggregate, ms_elementary, pias, na.rm = TRUE, r = -1)
++   Map(aggregate, ms_elementary, pias, na.rm = TRUE, order = -1)
 + )
 > 
 > paasche
diff --git a/vignettes/index-number-formulas.qmd b/vignettes/index-number-formulas.qmd
index 955ae15..639687d 100644
--- a/vignettes/index-number-formulas.qmd
+++ b/vignettes/index-number-formulas.qmd
@@ -52,7 +52,7 @@ corresponds to an arithmetic mean.
 
 ```{r}
 prices |>
-  elementary_index(price / back_price ~ period + business, r = 1)
+  elementary_index(price / back_price ~ period + business, order = 1)
 ```
 
 The Coggeshall index (equally-weighted harmonic mean of price relatives)
@@ -62,7 +62,7 @@ an order `r` of -1.
 
 ```{r}
 prices |>
-  elementary_index(price / back_price ~ period + business, r = -1)
+  elementary_index(price / back_price ~ period + business, order = -1)
 ```
 
 Weights can be added to make, for example, elementary indexes using the
diff --git a/vignettes/superlative-aggregation.qmd b/vignettes/superlative-aggregation.qmd
index c56e732..1438e9f 100644
--- a/vignettes/superlative-aggregation.qmd
+++ b/vignettes/superlative-aggregation.qmd
@@ -92,7 +92,7 @@ paasche <- Map(
   paasche_pias,
   na.rm = TRUE,
   include_ea = FALSE,
-  r = -1
+  order = -1
 ) |>
   Reduce(stack, x = _)