diff --git a/DESCRIPTION b/DESCRIPTION
index 70906aa..e831331 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: effectplots
 Title: Effect Plots
-Version: 0.2.2
+Version: 0.2.3
 Authors@R: 
     person("Michael", "Mayer", , "mayermichael79@gmail.com", role = c("aut", "cre"))
 Description: High-performance implementation of various effect plots
diff --git a/NEWS.md b/NEWS.md
index aab2719..f45c067 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# effectplots 0.2.3
+
+- More unit tests.
+- Better code explanations for .ale().
+
 # effectplots 0.2.2
 
 ### Minor improvement
diff --git a/R/ale.R b/R/ale.R
index 43bbd6d..5c3b4a0 100644
--- a/R/ale.R
+++ b/R/ale.R
@@ -1,7 +1,7 @@
 #' Accumulated Local Effects (ALE)
 #'
 #' @description
-#' Calculates ALE for one or multiple continuous features specified by `X`.
+#' Calculates uncentered ALE for one or multiple continuous features specified by `X`.
 #'
 #' The concept of ALE was introduced in Apley et al. (2020) as an alternative to
 #' partial dependence (PD). The Ceteris Paribus clause behind PD is a blessing and
@@ -17,7 +17,11 @@
 #' observations falling into this bin. This is repeated for all bins,
 #' and the values are *accumulated*.
 #'
-#' ALE values are plotted against right bin breaks.
+#' This implementation closely follows the implementation of Apley et al. (2020).
+#'  Notably, we also plot the values at the bin breaks, not at the bin means.
+#'  The main difference to Apley is that we use uniform binning, not quantile binning.
+#'  For large bins, we sample `ale_bin_size` observations, a step that is not necessary
+#'  with quantile binning. Furthermore, we don't center the values to mean 0.
 #'
 #' @details
 #' The function is a convenience wrapper around [feature_effects()], which calls
@@ -43,7 +47,7 @@
 #' M |> plot()
 #'
 #' M2 <- ale(fit, v = colnames(iris)[-1], data = iris, breaks = 5)
-#' plot(M2, share_y = "all")  # Only continuous variables shown
+#' plot(M2, share_y = "all") # Only continuous variables shown
 ale <- function(object, ...) {
   UseMethod("ale")
 }
@@ -65,8 +69,7 @@ ale.default <- function(
     ale_n = 50000L,
     ale_bin_size = 200L,
     seed = NULL,
-    ...
-) {
+    ...) {
   feature_effects.default(
     object = object,
     v = v,
@@ -106,8 +109,7 @@ ale.ranger <- function(
     ale_n = 50000L,
     ale_bin_size = 200L,
     seed = NULL,
-    ...
-) {
+    ...) {
   if (is.null(pred_fun)) {
     pred_fun <- function(model, newdata, ...) {
       stats::predict(model, newdata, ...)$predictions
@@ -149,8 +151,7 @@ ale.explainer <- function(
     ale_n = 50000L,
     ale_bin_size = 200L,
     seed = NULL,
-    ...
-) {
+    ...) {
   ale.default(
     object = object[["model"]],
     v = v,
@@ -187,8 +188,7 @@ ale.H2OModel <- function(
     ale_n = 50000L,
     ale_bin_size = 200L,
     seed = NULL,
-    ...
-) {
+    ...) {
   if (!requireNamespace("h2o", quietly = TRUE)) {
     stop("Package 'h2o' not installed")
   }
@@ -227,9 +227,8 @@ ale.H2OModel <- function(
 #' Per bin, the local effect \eqn{D_j} is calculated, and then accumulated over bins.
 #' \eqn{D_j} equals the difference between the partial dependence at the
 #' lower and upper bin breaks using only observations within bin.
-#' To plot the values, we can make a line plot of the resulting vector against
-#' upper bin breaks. Alternatively, the vector can be extended
-#' from the left by the value 0, and then plotted against *all* breaks.
+#' The values are to be plotted against bin breaks.
+#' Note that no centering is applied, i.e., the first value starts at 0.
 #'
 #' @param v Variable name in `data` to calculate ALE.
 #' @param data Matrix or data.frame.
@@ -242,7 +241,7 @@ ale.H2OModel <- function(
 #' @param g For internal use. The result of `as.factor(findInterval(...))`.
 #'   By default `NULL`.
 #' @inheritParams feature_effects
-#' @returns Vector representing one ALE per bin.
+#' @returns Vector representing one value per break.
 #' @export
 #' @seealso [partial_dependence()]
 #' @inherit ale references
@@ -262,47 +261,52 @@ ale.H2OModel <- function(
     bin_size = 200L,
     w = NULL,
     g = NULL,
-    ...
-) {
+    ...) {
   if (is.null(g)) {
     x <- if (is.data.frame(data)) data[[v]] else data[, v]
     g <- findInterval(
-      x, vec = breaks, rightmost.closed = TRUE, left.open = right, all.inside = TRUE
+      x = x, vec = breaks, rightmost.closed = TRUE, left.open = right, all.inside = TRUE
     )
     g <- collapse::qF(g, sort = FALSE)
   }
 
-  # List of bin indices. We remove empty or NA bins.
+  # List containing selected row indices per unsorted(!) bin ID
   J <- lapply(
     collapse::gsplit(g = g, use.g.names = TRUE),
     function(z) if (length(z) <= bin_size) z else sample(z, size = bin_size)
   )
+  # Remove empty or NA bins
   ok <- !is.na(names(J)) & lengths(J, use.names = FALSE) > 0L
   if (!all(ok)) {
     J <- J[ok]
   }
 
   # Before flattening the list J, we store bin counts
-  bin_n <- lengths(J, use.names = FALSE)
-  ix <- as.integer(names(J))
+  bin_sizes <- lengths(J, use.names = FALSE)
+  ix <- as.integer(names(J)) # Unsorted bin IDs
   J <- unlist(J, recursive = FALSE, use.names = FALSE)
 
-  # Empty bins will get an incremental effect of 0
-  out <- numeric(length(breaks) - 1L)
+  # Initialize local effects with 0
+  out <- numeric(length(breaks)) # first value corresponds to first left break
 
-  # Now we create a single prediction dataset. Lower bin edges first, then upper ones.
+  # Create single prediction data set. Lower bin edges first, then upper ones
+  # This makes the code harder to read, but is more efficient
   data_long <- collapse::ss(data, rep.int(J, 2L))
-  grid_long <- rep.int(c(breaks[ix], breaks[ix + 1L]), times = c(bin_n, bin_n))
+  grid_long <- rep.int(c(breaks[ix], breaks[ix + 1L]), times = c(bin_sizes, bin_sizes))
   if (is.data.frame(data_long)) {
     data_long[[v]] <- grid_long
   } else {
     data_long[, v] <- grid_long
   }
+
   pred <- prep_pred(
-    pred_fun(object, data_long, ...), trafo = trafo, which_pred = which_pred
+    pred_fun(object, data_long, ...),
+    trafo = trafo, which_pred = which_pred
   )
+
+  # Aggregate individual local effects
   n <- length(J)
-  out[ix] <- collapse::fmean(
+  out[ix + 1L] <- collapse::fmean(
     pred[(n + 1L):(2L * n)] - pred[1L:n],
     g = collapse::fdroplevels(g[J]),
     w = if (!is.null(w)) w[J],
diff --git a/R/feature_effects.R b/R/feature_effects.R
index 8bb8c6d..7f0bb55 100644
--- a/R/feature_effects.R
+++ b/R/feature_effects.R
@@ -594,7 +594,7 @@ calculate_stats <- function(
         w = ale_data$w,
         g = if (is.null(ale_data$ix)) ix else ix[ale_data$ix],
         ...
-      )
+      )[-1L] # drop value at first break as we have one value too much
       ok <- !is.na(out$bin_mid)
 
       # Centering possible?
diff --git a/man/ale.Rd b/man/ale.Rd
index 4050bf0..78b5f73 100644
--- a/man/ale.Rd
+++ b/man/ale.Rd
@@ -159,7 +159,7 @@ feature is discrete or continuous. This attribute might be lost when you manuall
 modify the data.frames.
 }
 \description{
-Calculates ALE for one or multiple continuous features specified by \code{X}.
+Calculates uncentered ALE for one or multiple continuous features specified by \code{X}.
 
 The concept of ALE was introduced in Apley et al. (2020) as an alternative to
 partial dependence (PD). The Ceteris Paribus clause behind PD is a blessing and
@@ -176,7 +176,11 @@ partial dependence difference between lower and upper bin break, using only
 observations falling into this bin. This is repeated for all bins,
 and the values are \emph{accumulated}.
 
-ALE values are plotted against right bin breaks.
+This implementation closely follows the implementation of Apley et al. (2020).
+Notably, we also plot the values at the bin breaks, not at the bin means.
+The main difference to Apley is that we use uniform binning, not quantile binning.
+For large bins, we sample \code{ale_bin_size} observations, a step that is not necessary
+with quantile binning. Furthermore, we don't center the values to mean 0.
 }
 \details{
 The function is a convenience wrapper around \code{\link[=feature_effects]{feature_effects()}}, which calls
@@ -199,7 +203,7 @@ M <- ale(fit, v = "Petal.Length", data = iris)
 M |> plot()
 
 M2 <- ale(fit, v = colnames(iris)[-1], data = iris, breaks = 5)
-plot(M2, share_y = "all")  # Only continuous variables shown
+plot(M2, share_y = "all") # Only continuous variables shown
 }
 \references{
 Apley, Daniel W., and Jingyu Zhu. 2020. \emph{Visualizing the Effects of Predictor Variables in Black Box Supervised Learning Models.}
diff --git a/man/dot-ale.Rd b/man/dot-ale.Rd
index 91d2598..c8f0981 100644
--- a/man/dot-ale.Rd
+++ b/man/dot-ale.Rd
@@ -55,16 +55,15 @@ By default \code{NULL}.}
 a \code{glm()} or (typically) \code{prob = TRUE} in classification models.}
 }
 \value{
-Vector representing one ALE per bin.
+Vector representing one value per break.
 }
 \description{
 This is a barebone implementation of Apley's ALE.
 Per bin, the local effect \eqn{D_j} is calculated, and then accumulated over bins.
 \eqn{D_j} equals the difference between the partial dependence at the
 lower and upper bin breaks using only observations within bin.
-To plot the values, we can make a line plot of the resulting vector against
-upper bin breaks. Alternatively, the vector can be extended
-from the left by the value 0, and then plotted against \emph{all} breaks.
+The values are to be plotted against bin breaks.
+Note that no centering is applied, i.e., the first value starts at 0.
 }
 \examples{
 fit <- lm(Sepal.Length ~ ., data = iris)
diff --git a/packaging.R b/packaging.R
index 53601b3..9284f7b 100644
--- a/packaging.R
+++ b/packaging.R
@@ -1,6 +1,6 @@
-#=============================================================================
+# =============================================================================
 # Put together the package
-#=============================================================================
+# =============================================================================
 
 # WORKFLOW: UPDATE EXISTING PACKAGE
 # 1) Modify package content and documentation.
@@ -15,7 +15,7 @@ library(usethis)
 use_description(
   fields = list(
     Title = "Effect Plots",
-    Version = "0.2.2",
+    Version = "0.2.3",
     Description = "High-performance implementation of various effect plots
     useful for regression and probabilistic classification tasks.
     The package includes partial dependence plots
@@ -47,7 +47,8 @@ use_gpl_license()
 
 # Your files that do not belong to the package itself (others are added by "use_* function")
 use_build_ignore(
-  c("^packaging.R$", "[.]Rproj$", "^logo.png$", "^claims.parquet$"), escape = FALSE
+  c("^packaging.R$", "[.]Rproj$", "^logo.png$", "^claims.parquet$"),
+  escape = FALSE
 )
 
 # Add short docu in Markdown (without running R code)
@@ -74,9 +75,9 @@ use_rcpp()
 # use_github_action("test-coverage")
 # use_github_action("pkgdown")
 
-#=============================================================================
+# =============================================================================
 # Finish package building (can use fresh session)
-#=============================================================================
+# =============================================================================
 
 library(devtools)
 
diff --git a/tests/testthat/test-ale.R b/tests/testthat/test-ale.R
index 38ed696..9050192 100644
--- a/tests/testthat/test-ale.R
+++ b/tests/testthat/test-ale.R
@@ -6,9 +6,68 @@ test_that("ale() is consistent with feature_effects()", {
   )
   suppressMessages(
     marg <- feature_effects(
-      fit, v = v, data = iris, calc_pred = FALSE, pd_n = 0, w = 1:150
+      fit,
+      v = v, data = iris, calc_pred = FALSE, pd_n = 0, w = 1:150
     )
   )
   expect_equal(ale, marg)
 })
 
+test_that(".ale() respects case weights", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+
+  res_w <- .ale(fit, v = v, data = iris, breaks = br, w = w)
+  res_uw <- .ale(fit, v = v, data = iris[ix, ], breaks = br)
+
+  expect_equal(res_w, res_uw)
+})
+
+test_that("ale() respects case weights", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+
+  res_w <- ale(fit, v = v, data = iris, breaks = br, w = w)[[1L]]
+  res_uw <- ale(fit, v = v, data = iris[ix, ], breaks = br)[[1L]]
+
+  expect_equal(res_w[-4L], res_uw[-4L])
+})
+
+test_that("the level order of g does not matter in .ale()", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  g1 <- factor(cut(iris[[v]], breaks = br, include.lowest = TRUE, labels = FALSE))
+  g2 <- factor(g1, levels = rev(levels(g1)))
+  res_1 <- .ale(fit, v = v, data = iris, breaks = br, g = g1)
+  res_2 <- .ale(fit, v = v, data = iris, breaks = br, g = g2)
+
+  expect_equal(res_1, res_2)
+})
+
+test_that(".ale() is consistent with uncentered ALEPlot v 1.1", {
+  # We use debugonce(ALEPlot) to see the uncentered values
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  K <- 5
+  br <- unique(stats::quantile(iris[[v]], probs = seq(0, 1, length.out = K + 1)))
+
+  # debugonce(ALEPlot)
+  # ALEPlot(
+  #   iris,
+  #   X.model = fit,
+  #   pred.fun = function(X.model, newdata) as.numeric(predict(X.model, newdata)),
+  #   J = v,
+  #   K = 5
+  # )
+  reference <- c(0.0000000, 0.6103576, 0.8673605, 0.9488453, 1.1848637, 1.9006788)
+  result <- .ale(fit, v = v, data = iris, breaks = br)
+
+  expect_equal(result, reference, tolerance = 1e-6)
+})
diff --git a/tests/testthat/test-average_observed.R b/tests/testthat/test-average_observed.R
index bb09f0a..a5b92c3 100644
--- a/tests/testthat/test-average_observed.R
+++ b/tests/testthat/test-average_observed.R
@@ -23,3 +23,25 @@ test_that("single vector input works", {
   expect_equal(out1, out2)
   expect_equal(out1, out3)
 })
+
+test_that("case weights are respected", {
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+  y <- iris$Sepal.Length
+
+  res_w <- average_observed(
+    iris[v],
+    y = y,
+    breaks = br,
+    w = w
+  )[[1L]]
+  res_uw <- average_observed(
+    iris[ix, v],
+    y = y[ix],
+    breaks = br
+  )[[1L]]
+
+  expect_equal(res_w[-4L], res_uw[-4L])
+})
diff --git a/tests/testthat/test-average_predicted.R b/tests/testthat/test-average_predicted.R
index 02f0a3b..d8616f9 100644
--- a/tests/testthat/test-average_predicted.R
+++ b/tests/testthat/test-average_predicted.R
@@ -16,3 +16,25 @@ test_that("single vector input works", {
   expect_equal(out1, out2)
   expect_equal(out1, out3)
 })
+
+test_that("case weights are respected", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+
+  res_w <- average_predicted(
+    iris[v],
+    pred = predict(fit, iris),
+    breaks = br,
+    w = w
+  )[[1L]]
+  res_uw <- average_predicted(
+    iris[ix, v],
+    pred = predict(fit, iris[ix, ]),
+    breaks = br
+  )[[1L]]
+
+  expect_equal(res_w[-4L], res_uw[-4L])
+})
diff --git a/tests/testthat/test-bias.R b/tests/testthat/test-bias.R
index 86eef94..1f2b4f7 100644
--- a/tests/testthat/test-bias.R
+++ b/tests/testthat/test-bias.R
@@ -4,7 +4,8 @@ v <- c("Sepal.Width", "Species")
 test_that("bias() is consistent with feature_effects()", {
   b <- bias(iris[v], resid = iris$Sepal.Length - predict(fit, iris))
   marg <- feature_effects(
-    fit, v = v, data = iris, y = iris$Sepal.Length, pd_n = 0, ale_n = 0
+    fit,
+    v = v, data = iris, y = iris$Sepal.Length, pd_n = 0, ale_n = 0
   )
   for (z in v) {
     out <- b[[z]]
@@ -22,3 +23,26 @@ test_that("single vector input works", {
   expect_equal(out1, out2)
   expect_equal(out1, out3)
 })
+
+test_that("case weights are respected", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+  r <- iris$Sepal.Length - predict(fit, iris)
+
+  res_w <- bias(
+    iris[v],
+    resid = r,
+    breaks = br,
+    w = w
+  )[[1L]]
+  res_uw <- bias(
+    iris[ix, v],
+    resid = r[ix],
+    breaks = br
+  )[[1L]]
+
+  expect_equal(res_w[-4L], res_uw[-4L])
+})
diff --git a/tests/testthat/test-partial_dependence.R b/tests/testthat/test-partial_dependence.R
index ad89e00..10b0690 100644
--- a/tests/testthat/test-partial_dependence.R
+++ b/tests/testthat/test-partial_dependence.R
@@ -3,7 +3,8 @@ test_that("partial_dependence() is consistent with feature_effects()", {
   v <- c("Sepal.Width", "Species")
   pd <- partial_dependence(fit, v = v, data = iris, w = 1:150)
   marg <- feature_effects(
-    fit, v = v, data = iris, calc_pred = FALSE, ale_n = 0, w = 1:150
+    fit,
+    v = v, data = iris, calc_pred = FALSE, ale_n = 0, w = 1:150
   )
   expect_equal(pd, marg)
 })
@@ -43,7 +44,34 @@ test_that(".pd() respects ... argument for predictions", {
 
   pd1 <- .pd(fit2, v = "Species", data = iris, grid = unique(iris$Species))
   pd2 <- .pd(
-    fit2, v = "Species", data = iris, grid = unique(iris$Species), type = "response"
+    fit2,
+    v = "Species", data = iris, grid = unique(iris$Species), type = "response"
   )
   expect_false(identical(pd1, pd2))
 })
+
+test_that(".pd() respects case weights", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+
+  res_w <- .pd(fit, v = v, data = iris, grid = br, w = w)
+  res_uw <- .pd(fit, v = v, data = iris[ix, ], grid = br)
+
+  expect_equal(res_w, res_uw)
+})
+
+test_that("partial_dependence() respects case weights", {
+  fit <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
+  v <- "Sepal.Width"
+  br <- c(2, 2.5, 3, 3.5, 4.5)
+  w <- c(rep(1L, times = 100L), rep(2L, times = 50L))
+  ix <- rep(1:nrow(iris), times = w)
+
+  res_w <- partial_dependence(fit, v = v, data = iris, breaks = br, w = w)[[1L]]
+  res_uw <- partial_dependence(fit, v = v, data = iris[ix, ], breaks = br)[[1L]]
+
+  expect_equal(res_w[-4L], res_uw[-4L])
+})