diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 1a2b762377..06cffc97d9 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,8 +1,6 @@
 name: CI
 on:
   pull_request:
-    branches:
-      - master
   push:
     tags:
       - '*'
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 54ef752bc9..eb5fb57b0e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## Unreleased
 
+- Added `aggregate` transformation for flexible data aggregation with automatic grouping, custom labels, and support for scalar and vector-valued aggregations [#696](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/696).
+
 ## v0.11.9 - 2025-10-10
 
 - Improved error message when two layers with incompatible continuous data are combined [#692](https://github.com/MakieOrg/AlgebraOfGraphics.jl/pull/692).
diff --git a/docs/src/reference/analyses.md b/docs/src/reference/analyses.md
index de3567eb39..7966fc01fc 100644
--- a/docs/src/reference/analyses.md
+++ b/docs/src/reference/analyses.md
@@ -10,35 +10,35 @@ EditURL = "analyses.jl"
 histogram
 ```
 
-````@example analyses
+```@example analyses
 using AlgebraOfGraphics, CairoMakie
 set_aog_theme!()
 
 df = (x=randn(5000), y=randn(5000), z=rand(["a", "b", "c"], 5000))
 specs = data(df) * mapping(:x, layout=:z) * histogram(bins=range(-2, 2, length=15))
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, dodge=:z, color=:z) * histogram(bins=range(-2, 2, length=15))
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, stack=:z, color=:z) * histogram(bins=range(-2, 2, length=15))
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) *
     mapping((:x, :z) => ((x, z) -> x + 5 * (z == "b")) => "new x", col=:z) *
     histogram(datalimits=extrema, bins=20)
 draw(specs, facet=(linkxaxes=:minimal,))
-````
+```
 
-````@example analyses
+```@example analyses
 data(df) * mapping(:x, :y, layout=:z) * histogram(bins=15) |> draw
-````
+```
 
 ## Density
 
@@ -46,38 +46,38 @@ data(df) * mapping(:x, :y, layout=:z) * histogram(bins=15) |> draw
 AlgebraOfGraphics.density
 ```
 
-````@example analyses
+```@example analyses
 df = (x=randn(5000) .+ repeat([0, 2, 4, 6], inner = 1250), y=randn(5000), z=repeat(["a", "b", "c", "d"], inner = 1250))
 specs = data(df) * mapping(:x, layout=:z) * AlgebraOfGraphics.density()
 
 draw(specs)
-````
+```
 
 ```@example analyses
 data(df) * mapping(:x, layout=:z) * AlgebraOfGraphics.density(datalimits = (0, 8)) |> draw
 ```
 
-````@example analyses
+```@example analyses
 draw(specs * visual(direction = :y))
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) *
     mapping((:x, :z) => ((x, z) -> x + 5 * (z ∈ ["b", "d"])) => "new x", layout=:z) *
     AlgebraOfGraphics.density(datalimits=extrema)
 draw(specs, facet=(linkxaxes=:minimal,))
-````
+```
 
-````@example analyses
+```@example analyses
 data(df) * mapping(:x, :y, layout=:z) * AlgebraOfGraphics.density(npoints=50) |> draw
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, :y, layout=:z) *
     AlgebraOfGraphics.density(npoints=50) * visual(Surface)
 
 draw(specs, axis=(type=Axis3, zticks=0:0.1:0.2, limits=(nothing, nothing, (0, 0.2))))
-````
+```
 
 ## Frequency
 
@@ -85,21 +85,21 @@ draw(specs, axis=(type=Axis3, zticks=0:0.1:0.2, limits=(nothing, nothing, (0, 0.
 frequency
 ```
 
-````@example analyses
+```@example analyses
 df = (x=rand(["a", "b", "c"], 100), y=rand(["a", "b", "c"], 100), z=rand(["a", "b", "c"], 100))
 specs = data(df) * mapping(:x, layout=:z) * frequency()
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, layout=:z, color=:y, stack=:y) * frequency()
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, :y, layout=:z) * frequency()
 draw(specs)
-````
+```
 
 ## Expectation
 
@@ -107,21 +107,21 @@ draw(specs)
 expectation
 ```
 
-````@example analyses
+```@example analyses
 df = (x=rand(["a", "b", "c"], 100), y=rand(["a", "b", "c"], 100), z=rand(100), c=rand(["a", "b", "c"], 100))
 specs = data(df) * mapping(:x, :z, layout=:c) * expectation()
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, :z, layout=:c, color=:y, dodge=:y) * expectation()
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 specs = data(df) * mapping(:x, :y, :z, layout=:c) * expectation()
 draw(specs)
-````
+```
 
 ## Linear
 
@@ -129,14 +129,14 @@ draw(specs)
 linear
 ```
 
-````@example analyses
+```@example analyses
 x = 1:0.05:10
 a = rand(1:7, length(x))
 y = 1.2 .* x .+ a .+ 0.5 .* randn.()
 df = (; x, y, a)
 specs = data(df) * mapping(:x, :y, color=:a => nonnumeric) * (linear() + visual(Scatter))
 draw(specs)
-````
+```
 
 ## Smoothing
 
@@ -144,14 +144,14 @@ draw(specs)
 smooth
 ```
 
-````@example analyses
+```@example analyses
 x = 1:0.05:10
 a = rand(1:7, length(x))
 y = sin.(x) .+ a .+ 0.1 .* randn.()
 df = (; x, y, a)
 specs = data(df) * mapping(:x, :y, color=:a => nonnumeric) * (smooth() + visual(Scatter))
 draw(specs)
-````
+```
 
 ## Contours
 
@@ -159,23 +159,23 @@ draw(specs)
 contours
 ```
 
-````@example analyses
+```@example analyses
 x = repeat(1:10, 10)
 y = repeat(11:20, inner = 10)
 z = sqrt.(x .* y)
 df = (; x, y, z)
 specs = data(df) * mapping(:x, :y, :z) * contours(levels = 8)
 draw(specs)
-````
+```
 
-````@example analyses
+```@example analyses
 x = repeat(1:10, 10)
 y = repeat(11:20, inner = 10)
 z = sqrt.(x .* y)
 df = (; x, y, z)
 specs = data(df) * mapping(:x, :y, :z) * contours(levels = 8, labels = true)
 draw(specs)
-````
+```
 
 ## Filled Contours
 
@@ -205,3 +205,93 @@ specs = data(df) *
     filled_contours(levels = [-Inf, 5, 8, 10, 12, 13, 14, Inf])
 draw(specs, scales(Color = (; palette = clipped(from_continuous(:plasma), low = :cyan, high = :red))))
 ```
+
+## Aggregate
+
+```@docs
+aggregate
+```
+
+The `aggregate` transformation allows you to perform flexible aggregations on your data.
+All mapped columns that are not explicitly aggregated are automatically used for grouping.
+
+This analysis layer is intended for aggregations that are only needed for a visualization, otherwise it may make more sense to compute values in a separate data wrangling step and add a separate `data` layer.
+
+### Basic Aggregation
+
+Compute the mean body mass for each penguin species:
+
+```@example analyses
+using AlgebraOfGraphics
+using Statistics
+
+penguins = AlgebraOfGraphics.penguins()
+
+data(penguins) *
+    mapping(:species, :body_mass_g) *
+    aggregate(2 => mean) *
+    visual(BarPlot) |> draw
+```
+
+### Multiple Grouping Dimensions
+
+Group by both species and sex, computing mean body mass:
+
+```@example analyses
+data(penguins) *
+    mapping(:species, :body_mass_g, color = :sex, dodge = :sex) *
+    aggregate(2 => mean) *
+    visual(BarPlot) |> draw
+```
+
+### Aggregating Multiple Columns
+
+Compute both mean and standard deviation:
+
+```@example analyses
+data(penguins) *
+    mapping(:species, :body_mass_g) *
+    (
+        aggregate(2 => mean) * visual(BarPlot) +
+        aggregate(2 => mean, 2 => std => 3) * visual(Errorbars)
+    ) |> draw
+```
+
+### Splitting Aggregation Results
+
+Sometimes an aggregation function may return multiple values which
+should form separate inputs for the subsequent visual. In this
+case you can assign a vector of accessor specifications via the pair syntax.
+Each vector element must specify an accessor function (here `first` and `last`) and the mapping that the result should be assigned to, either given as integers for positional arguments or symbols for named arguments.
+
+```@example analyses
+data(penguins) *
+    mapping(:species, :body_mass_g, color = :sex, dodge_x = :sex) *
+    aggregate(2 => extrema => [first => 2, last => 3]) *
+    visual(Rangebars, linewidth = 3) |> draw(scales(DodgeX = (; width = 0.2)))
+```
+
+### Vector-Valued Aggregations
+
+Aggregate functions can return vectors, which will be automatically expanded. If you have another aggregation that returns scalars, the scalars will be repeated to match the length of the vector aggregation. If you have multiple aggregations returning vectors, the lengths of all vectors in a given group must be the same.
+
+```@example analyses
+# Get lower and upper quartiles as a vector
+lower_upper_quartile(x) = quantile(x, [0.25, 0.75])
+
+data(penguins) *
+    mapping(:species, :body_mass_g, color = :sex) *
+    aggregate(2 => lower_upper_quartile) *
+    visual(Scatter, markersize = 15) |> draw
+```
+
+### Custom Labels
+
+Provide custom labels for aggregated outputs:
+
+```@example analyses
+data(penguins) *
+    mapping(:species, :body_mass_g) *
+    aggregate(2 => mean => "Average Mass (g)") *
+    visual(BarPlot) |> draw
+```
diff --git a/src/AlgebraOfGraphics.jl b/src/AlgebraOfGraphics.jl
index 93774a15de..aee563b9d9 100644
--- a/src/AlgebraOfGraphics.jl
+++ b/src/AlgebraOfGraphics.jl
@@ -35,7 +35,7 @@ export hideinnerdecorations!, deleteemptyaxes!
 export Layer, Layers, ProcessedLayer, ProcessedLayers, zerolayer
 export Entry, AxisEntries
 export renamer, sorter, nonnumeric, verbatim, presorted
-export density, histogram, linear, smooth, expectation, frequency, contours, filled_contours
+export density, histogram, linear, smooth, expectation, frequency, contours, filled_contours, aggregate, highlight
 export visual, data, geodata, dims, mapping
 export datetimeticks
 export draw, draw!
@@ -72,8 +72,10 @@ include("transformations/histogram.jl")
 include("transformations/groupreduce.jl")
 include("transformations/frequency.jl")
 include("transformations/expectation.jl")
+include("transformations/aggregate.jl")
 include("transformations/contours.jl")
 include("transformations/filled_contours.jl")
+include("transformations/highlight.jl")
 include("guides/guides.jl")
 include("guides/legend.jl")
 include("guides/colorbar.jl")
diff --git a/src/algebra/layer.jl b/src/algebra/layer.jl
index d604f48d64..c35764d7be 100644
--- a/src/algebra/layer.jl
+++ b/src/algebra/layer.jl
@@ -169,8 +169,9 @@ end
 
 unnest(vs::AbstractArray, indices) = map(k -> [el[k] for el in vs], indices)
 
-unnest_arrays(vs) = unnest(vs, keys(first(vs)))
+unnest_arrays(vs) = isempty(vs) ? [[]] : unnest(vs, keys(first(vs)))
 function unnest_dictionaries(vs)
+    isempty(vs) && return Dictionary()
     return Dictionary(Dict((k => [el[k] for el in vs] for k in collect(keys(first(vs))))))
 end
 slice(v, c) = map(el -> getnewindex(el, c), v)
@@ -192,6 +193,21 @@ function Base.map(f, processedlayer::ProcessedLayer)
     return ProcessedLayer(processedlayer; positional, named)
 end
 
+function filtermap(f, processedlayer::ProcessedLayer)
+    axs = shape(processedlayer)
+    outputs = map(CartesianIndices(axs)) do c
+        return f(slice(processedlayer.positional, c), slice(processedlayer.named, c))
+    end
+    to_remove = outputs .== nothing
+    deleteat!(outputs, to_remove)
+    primary = map(processedlayer.primary) do value
+        value[.!to_remove]
+    end
+    
+    positional, named = unnest_arrays(map(first, outputs)), unnest_dictionaries(map(last, outputs))
+    return ProcessedLayer(processedlayer; positional, named, primary)
+end
+
 ## Get scales from a `ProcessedLayer`
 
 function uniquevalues(v::AbstractArray)
@@ -200,7 +216,7 @@ function uniquevalues(v::AbstractArray)
     return collect(uniquesorted(_v, perm))
 end
 
-to_label(label::AbstractString) = label
+to_label(label) = label
 to_label(labels::AbstractArray) = reduce(mergelabels, labels)
 
 # merge dict2 into dict but translate keys first using remapdict
diff --git a/src/guides/legend.jl b/src/guides/legend.jl
index 5c32be3a3b..f971f77a36 100644
--- a/src/guides/legend.jl
+++ b/src/guides/legend.jl
@@ -96,6 +96,9 @@ end
 
 categorical_scales_mergeable(c1, c2) = false # there can be continuous scales in the mix, like markersize
 
+is_empty_categorical_scale(s::CategoricalScale) = isempty(datavalues(s))
+is_empty_categorical_scale(s::ContinuousScale) = false
+
 function compute_legend(grid::Matrix{<:Union{AxisEntries, AxisSpecEntries}}; order::Union{Nothing, AbstractVector})
     # gather valid named scales
     scales_categorical = legendable_scales(Val(:categorical), first(grid).categoricalscales)
@@ -103,6 +106,9 @@ function compute_legend(grid::Matrix{<:Union{AxisEntries, AxisSpecEntries}}; ord
 
     scales = Iterators.flatten((pairs(scales_categorical), pairs(scales_continuous)))
 
+    # if no legendable scale is present, return nothing
+    isempty(scales) || all(x -> all(is_empty_categorical_scale, last(x)), scales) && return nothing
+
     processedlayers = first(grid).processedlayers
 
     # we can't loop over all processedlayers here because one layer can be sliced into multiple processedlayers
diff --git a/src/scales.jl b/src/scales.jl
index 3ee3b2a31e..ed9b61163d 100644
--- a/src/scales.jl
+++ b/src/scales.jl
@@ -415,7 +415,7 @@ continuous_aes_props(type::Type{<:Aesthetic}, props_dict::Dictionary{Symbol, Any
 
 struct ContinuousScale{T}
     extrema::NTuple{2, T}
-    label::Union{AbstractString, Nothing}
+    label # nothing or any type workable as a label
     force::Bool
     props::ContinuousScaleProps
 end
diff --git a/src/transformations/aggregate.jl b/src/transformations/aggregate.jl
new file mode 100644
index 0000000000..8559e62df8
--- /dev/null
+++ b/src/transformations/aggregate.jl
@@ -0,0 +1,508 @@
+struct AggregationOutput
+    accessor::Union{Nothing, Base.Callable}  # Function to extract from aggregation result (nothing means use as-is)
+    destination::Union{Int, Symbol}  # Where to place this output
+    label::Any  # Optional label (can be String, RichText, etc.)
+    scaleid::Union{Nothing, ScaleID}  # Optional scale id
+end
+
+struct ParsedAggregation
+    target::Union{Int, Symbol}  # Source column to aggregate
+    aggfunc::Base.Callable  # Aggregation function
+    outputs::Vector{AggregationOutput}  # One or more outputs from this aggregation
+end
+
+struct AggregateAnalysis{A, G}
+    aggregations::A
+    groupby::G
+end
+
+# Parse a single output spec: destination or destination => label or destination => label => scale_id
+function _parse_output_spec(dest::Union{Int, Symbol})
+    return AggregationOutput(nothing, dest, nothing, nothing)
+end
+
+# Helper to parse the "rest" part after destination
+_parse_label_and_scale(label) = (label, nothing)
+_parse_label_and_scale(p::Pair{<:Any, ScaleID}) = (first(p), last(p))
+
+# dest => something (dispatch on the something)
+function _parse_output_spec(p::Pair)
+    dest = first(p)
+    rest = last(p)
+    label, scaleid = _parse_label_and_scale(rest)
+    return AggregationOutput(nothing, dest, label, scaleid)
+end
+
+# Parse split spec: accessor => output_spec
+function _parse_split(p::Pair{<:Base.Callable, <:Union{Int, Symbol}})
+    accessor = first(p)
+    dest = last(p)
+    return AggregationOutput(accessor, dest, nothing, nothing)
+end
+
+function _parse_split(p::Pair{<:Base.Callable, <:Pair})
+    accessor = first(p)
+    dest_spec = last(p)
+    output = _parse_output_spec(dest_spec)
+    return AggregationOutput(accessor, output.destination, output.label, output.scaleid)
+end
+
+# Fallback for when the array element type is inferred as Pair{Function, Any}
+# This delegates to internal dispatch methods
+function _parse_split(p::Pair{<:Union{Function, Type}})
+    accessor = first(p)
+    dest_spec = last(p)
+    return _parse_split_internal(accessor, dest_spec)
+end
+
+# Internal dispatch for split parsing
+_parse_split_internal(accessor, dest::Union{Int, Symbol}) =
+    AggregationOutput(accessor, dest, nothing, nothing)
+
+function _parse_split_internal(accessor, dest_spec::Pair)
+    output = _parse_output_spec(dest_spec)
+    return AggregationOutput(accessor, output.destination, output.label, output.scaleid)
+end
+
+# Helper to parse aggregation spec: just a function
+_parse_agg_spec(target, f::Base.Callable) =
+    ParsedAggregation(target, f, [AggregationOutput(nothing, target, nothing, nothing)])
+
+# Helper to parse the splits/label part of function => splits_or_label
+_parse_outputs(target, aggfunc, splits::AbstractVector) =
+    ParsedAggregation(target, aggfunc, map(_parse_split, splits))
+
+# When the third element is an Int or Symbol, treat it as a destination
+function _parse_outputs(target, aggfunc, dest::Union{Int, Symbol})
+    return ParsedAggregation(target, aggfunc, [AggregationOutput(nothing, dest, nothing, nothing)])
+end
+
+# Otherwise it's a label and/or scale_id
+function _parse_outputs(target, aggfunc, label_and_or_scale)
+    label, scaleid = _parse_label_and_scale(label_and_or_scale)
+    return ParsedAggregation(target, aggfunc, [AggregationOutput(nothing, target, label, scaleid)])
+end
+
+# Helper to parse aggregation spec: function => something
+function _parse_agg_spec(target, p::Pair{<:Base.Callable})
+    aggfunc = first(p)
+    outputs_spec = last(p)
+    return _parse_outputs(target, aggfunc, outputs_spec)
+end
+
+"""
+    aggregate(aggregations...; named_aggregations...)
+
+Perform flexible aggregation of data. Specify which columns to aggregate explicitly;
+all other mapped columns are automatically used for grouping.
+
+# Arguments
+- Positional arguments are aggregation specifications in the form:
+  - `target => aggfunc` where target is an Int (positional) or Symbol (named)
+  - `target => aggfunc => dest` to place output at a different position
+  - `target => aggfunc => [accessor => dest, ...]` to split aggregation results
+- Named arguments are aggregation functions for named mappings (e.g., `color = mean`)
+
+# Labeling and Custom Scales
+You can customize labels and assign outputs to custom scales:
+- `target => aggfunc => label` - Set a custom label (if label is not an Int/Symbol)
+- `target => aggfunc => label => scale(:scaleid)` - Set label and assign to a custom scale
+- For split outputs: `accessor => dest => label` or `accessor => dest => label => scale(:scaleid)`
+
+# Examples
+
+```julia
+# Aggregate y values (position 2), group by x (position 1)
+data(...) * mapping(:time, :value) * 
+    aggregate(2 => median)
+
+# Aggregate and place in different position
+data(...) * mapping(:time, :value) *
+    aggregate(2 => mean, 2 => std => 3) *
+    visual(Errorbars)
+
+# Aggregate x values (position 1), group by y (position 2)
+data(...) * mapping(:value, :time) * 
+    aggregate(1 => mean)
+
+# Multiple aggregations with named argument
+data(...) * mapping(:time, :value, color=:group) * 
+    aggregate(2 => mean, color = length)
+
+# Group by multiple dimensions (x and y), aggregate z
+data(...) * mapping(:x, :y, :z) * 
+    aggregate(3 => mean)
+
+# Split extrema into separate positions for range bars
+data(...) * mapping(:x, :y) * 
+    aggregate(2 => extrema => [first => 2, last => 3]) *
+    visual(Rangebars)
+
+# Aggregate multiple columns
+data(...) * mapping(:x, :y1, :y2) *
+    aggregate(2 => mean, 3 => median)
+
+# Custom labels
+data(...) * mapping(:x, :y) *
+    aggregate(2 => mean => "Average Y")
+
+# Custom labels with LaTeX
+data(...) * mapping(:x, :y) *
+    aggregate(2 => mean => L"\\bar{y}")
+
+# Split outputs with custom labels and scale
+data(...) * mapping(:x, :y) *
+    aggregate(
+        2 => extrema => [
+            first => 2 => "Minimum",
+            last => :color => "Maximum" => scale(:color2)
+        ]
+    ) *
+    visual(Scatter) |>
+    draw(scales(color2 = (; colormap = :thermal)))
+
+# Custom scale for aggregated output
+data(...) * mapping(:x, :y, :z) *
+    aggregate(3 => sum => "Total" => scale(:mycolor)) *
+    visual(Heatmap) |>
+    draw(scales(mycolor = (; colormap = :viridis)))
+```
+"""
+function aggregate(args...; named_aggs...)
+    # All arguments should be aggregation specifications (target => aggfunc)
+    aggregations = Pair[]
+
+    for arg in args
+        # arg should be a Pair with target => aggfunc
+        if !(arg isa Pair)
+            throw(ArgumentError("Each positional argument to aggregate must be a Pair like `target => aggfunc`, got $(typeof(arg))"))
+        end
+        push!(aggregations, arg)
+    end
+
+    # Add named aggregations
+    for (name, aggfunc) in pairs(named_aggs)
+        push!(aggregations, name => aggfunc)
+    end
+
+    # No explicit groupby needed - it will be inferred from what's not aggregated
+    return transformation(AggregateAnalysis(Tuple(aggregations), nothing))
+end
+
+function (a::AggregateAnalysis)(input::ProcessedLayer)
+    N = length(input.positional)
+
+    # Collect all targets being aggregated
+    aggregated_targets = Set{Union{Int, Symbol}}()
+    for (target, _) in a.aggregations
+        push!(aggregated_targets, target)
+    end
+
+    # Infer grouping columns: all positional indices not being aggregated
+    grouping_indices = Int[]
+    for i in 1:N
+        if !(i in aggregated_targets)
+            push!(grouping_indices, i)
+        end
+    end
+
+    # Also auto-group by named arguments that aren't being aggregated
+    grouping_names = Symbol[]
+    for key in keys(input.named)
+        if !(key in aggregated_targets)
+            push!(grouping_names, key)
+        end
+    end
+
+    # Build summaries for unique values in grouping dimensions
+    summaries = [
+        mapreduce(collect ∘ uniquesorted, mergesorted, input.positional[idx])
+            for idx in grouping_indices
+    ]
+
+    # Parse all aggregation specs once and compute labels for each output
+    parsed_aggregations = map(a.aggregations) do (target, agg_spec)
+        parsed = _parse_agg_spec(target, agg_spec)
+
+        # Get original label for this target
+        original_label = get(input.labels, target, "")
+
+        # Generate labels for each output
+        func_name = string(nameof(parsed.aggfunc))
+        base_label = isempty(original_label) ? func_name : "$(func_name)($(original_label))"
+
+        # Update outputs with generated labels where not provided
+        outputs = map(parsed.outputs) do output
+            if output.label !== nothing
+                # User provided explicit label
+                return output
+            elseif output.accessor !== nothing
+                # Generate label with accessor name
+                accessor_name = string(nameof(output.accessor))
+                generated_label = "$(accessor_name)($(base_label))"
+                return AggregationOutput(output.accessor, output.destination, generated_label, output.scaleid)
+            else
+                # Use base label
+                return AggregationOutput(output.accessor, output.destination, base_label, output.scaleid)
+            end
+        end
+
+        return ParsedAggregation(target, parsed.aggfunc, outputs)
+    end
+
+    # Build output labels dictionary (Any to support RichText, String, etc.)
+    # Wrap labels in fill() to make them broadcastable
+    # Also build scale_mapping dictionary to map positions/names to custom scale ids
+    output_labels = Dict{Union{Int, Symbol}, Any}()
+    scale_mapping = Dictionary{KeyType, Symbol}()
+
+    for parsed in parsed_aggregations
+        for output in parsed.outputs
+            output_labels[output.destination] = fill(output.label)
+            if output.scaleid !== nothing
+                insert!(scale_mapping, output.destination, output.scaleid.id)
+            end
+        end
+    end
+
+    # Perform aggregations and build output in a single map over input
+    output = map(input) do p, n
+        # Extract grouping keys once (same for all aggregations)
+        # Include both positional and named grouping columns
+        positional_keys = Tuple(p[idx] for idx in grouping_indices)
+        named_keys = Tuple(n[name] for name in grouping_names)
+        grouping_key_columns = (positional_keys..., named_keys...)
+
+        # Handle case where there are no grouping columns (single group)
+        perm, group_perm, actual_keys = if isempty(grouping_key_columns)
+            # No grouping - treat all data as a single group
+            # Get number of rows from first positional column
+            n_rows = length(first(p))
+            # Create identity permutation and a single group containing all indices
+            perm = 1:n_rows
+            # GroupPerm expects a vector of index ranges for each group
+            # For a single group with all rows, that's just [1:n_rows]
+            group_ranges = [1:n_rows]
+            actual_keys = [()]  # Single empty key tuple
+            (perm, group_ranges, actual_keys)
+        else
+            sa = StructArray(map(fast_hashed, grouping_key_columns))
+            perm = sortperm(sa)
+            group_perm = GroupPerm(sa, perm)
+
+            # Extract actual keys that exist in the data (unhashed)
+            actual_keys = map(group_perm) do idxs
+                idx = perm[first(idxs)]
+                # Extract the unhashed key values for this group
+                return map(k -> k[idx], grouping_key_columns)
+            end
+            (perm, group_perm, actual_keys)
+        end
+
+        # Build output dictionary - will contain all results indexed by position or symbol
+        outputs = Dict{Union{Int, Symbol}, Any}()
+
+        # Process all aggregations first to determine if we need to expand groups
+        aggregation_results = Dict{Union{Int, Symbol}, Any}()
+        targets = Dict{Union{Int, Symbol}, Union{Int, Symbol}}()
+
+        for parsed in parsed_aggregations
+            target = parsed.target
+            aggfunc = parsed.aggfunc
+
+            # Validate target and extract target values
+            if target isa Integer
+                if target < 1 || target > N
+                    throw(ArgumentError("aggregation target $target out of bounds for $N positional arguments"))
+                end
+                if target in grouping_indices
+                    throw(ArgumentError("cannot aggregate positional argument $target which is used for grouping"))
+                end
+                target_values = p[target]
+            elseif target isa Symbol
+                if !haskey(n, target)
+                    throw(ArgumentError("aggregation target :$target not found in named arguments"))
+                end
+                if target in grouping_names
+                    throw(ArgumentError("cannot aggregate named argument :$target which is used for grouping"))
+                end
+                target_values = n[target]
+            else
+                throw(ArgumentError("aggregation target must be an Integer or Symbol, got $(typeof(target))"))
+            end
+
+            # Apply aggregation using the precomputed GroupPerm
+            result = map(group_perm) do idxs
+                group_indices = perm[idxs]
+                group_values = view(target_values, group_indices)
+                return aggfunc(group_values)
+            end
+
+            # Flatten multidimensional results
+            result = result isa AbstractArray && ndims(result) > 1 ? vec(result) : result
+
+            # Process each output from this aggregation
+            for output in parsed.outputs
+                destination = output.destination
+
+                if haskey(aggregation_results, destination)
+                    throw(ArgumentError("Output slot $(repr(destination)) of `aggregate` was assigned multiple times. By default, aggregated mappings are routed to their original position, for example `2 => std` is routed to positional arg 2. You can use the pair syntax `2 => std => 3` to route to a different positional or named argument."))
+                end
+
+                # Apply accessor if present, otherwise use result as-is
+                final_result = if output.accessor !== nothing
+                    map(output.accessor, result)
+                else
+                    result
+                end
+
+                aggregation_results[destination] = final_result
+                targets[destination] = target
+            end
+        end
+
+        # Detect which results are vector-valued and determine group lengths
+        # Check the element type of the result vectors
+        vector_valued_results = Dict{Union{Int, Symbol}, Bool}()
+        for (destination, result) in aggregation_results
+            # Check if the element type is a subtype of AbstractVector
+            element_type = eltype(result)
+            if element_type <: AbstractArray && !(element_type <: AbstractVector)
+                target = targets[destination]
+                dims = size(first(result))
+
+                # Build descriptive error message with column description
+                target_desc = if target isa Integer
+                    "positional argument $target"
+                else
+                    "argument :$target"
+                end
+
+                throw(ArgumentError("Aggregation of $target_desc returned $(length(dims))-dimensional arrays with size $dims. Only scalars or 1-dimensional vectors are supported."))
+            end
+            vector_valued_results[destination] = element_type <: AbstractVector
+        end
+
+        # Check if we have any vector results
+        has_vector_results = any(values(vector_valued_results))
+
+        if has_vector_results
+            # Vector-valued aggregation: determine lengths and validate consistency
+            # Compute the length of each group's result vector from vector-valued results
+            group_lengths = map(enumerate(actual_keys)) do (group_idx, key)
+                lengths_this_group = Int[]
+                for (destination, result) in aggregation_results
+                    if vector_valued_results[destination] && !isempty(result)
+                        push!(lengths_this_group, length(result[group_idx]))
+                    end
+                end
+
+                # Validate that all vector results for this group have the same length
+                if !isempty(lengths_this_group) && !allequal(lengths_this_group)
+                    vector_dests_lengths = [
+                        (d, length(aggregation_results[d][group_idx]))
+                            for (d, is_vec) in vector_valued_results
+                            if is_vec
+                    ]
+                    throw(ArgumentError("Inconsistent vector lengths for group at index $group_idx (key=$(key)): $vector_dests_lengths. All vector-valued aggregations must return the same length for each group."))
+                end
+
+                return isempty(lengths_this_group) ? 0 : first(lengths_this_group)
+            end
+
+            # Expand grouping columns by repeating keys according to their result lengths
+            # Handle positional grouping columns
+            for (i, group_idx) in enumerate(grouping_indices)
+                expanded = mapreduce(vcat, enumerate(actual_keys)) do (gi, key)
+                    fill(key[i], group_lengths[gi])
+                end
+                outputs[group_idx] = expanded
+            end
+
+            # Handle named grouping columns
+            num_positional_groups = length(grouping_indices)
+            for (i, group_name) in enumerate(grouping_names)
+                key_idx = num_positional_groups + i
+                expanded = mapreduce(vcat, enumerate(actual_keys)) do (gi, key)
+                    fill(key[key_idx], group_lengths[gi])
+                end
+                outputs[group_name] = expanded
+            end
+
+            # Process all aggregation results: concatenate vectors, expand scalars
+            for (destination, result) in aggregation_results
+                if vector_valued_results[destination]
+                    # Vector result: concatenate
+                    concatenated = mapreduce(vcat, result) do val
+                        val
+                    end
+                    outputs[destination] = concatenated
+                else
+                    # Scalar result: expand to match group lengths
+                    total_length = sum(group_lengths)
+                    expanded = similar(result, total_length)
+                    offset = 1
+                    for (gi, val) in enumerate(result)
+                        len = group_lengths[gi]
+                        fill!(view(expanded, offset:(offset + len - 1)), val)
+                        offset += len
+                    end
+                    outputs[destination] = expanded
+                end
+            end
+        else
+            # Scalar aggregation: use keys directly
+            # Handle positional grouping columns
+            for (i, group_idx) in enumerate(grouping_indices)
+                outputs[group_idx] = [key[i] for key in actual_keys]
+            end
+
+            # Handle named grouping columns
+            num_positional_groups = length(grouping_indices)
+            for (i, group_name) in enumerate(grouping_names)
+                key_idx = num_positional_groups + i
+                outputs[group_name] = [key[key_idx] for key in actual_keys]
+            end
+
+            # Use aggregation results as-is
+            for (destination, result) in aggregation_results
+                outputs[destination] = result
+            end
+        end
+
+        # Separate positional and named results
+        positional_keys = sort([k for k in keys(outputs) if k isa Integer])
+        named_keys = [k for k in keys(outputs) if k isa Symbol]
+
+        # Validate positional keys form a contiguous range from 1 to max
+        if !isempty(positional_keys)
+            max_pos = maximum(positional_keys)
+            expected = 1:max_pos
+            missing_positions = setdiff(expected, positional_keys)
+            if !isempty(missing_positions)
+                throw(ArgumentError("positional outputs must be contiguous, got $positional_keys missing $missing_positions"))
+            end
+        end
+
+        # Build positional array from sorted keys
+        positional = [outputs[k] for k in positional_keys]
+
+        # Build named arguments from symbol-keyed results
+        named_dict = Dictionary{Symbol, Any}()
+        for k in named_keys
+            insert!(named_dict, k, outputs[k])
+        end
+
+        named = merge(n, named_dict)
+
+        return positional, named
+    end
+
+    # Apply labels to the output
+    labels = set(output.labels, (k => v for (k, v) in output_labels)...)
+
+    # Merge scale_mapping with existing scale_mapping from output
+    merged_scale_mapping = merge(output.scale_mapping, scale_mapping)
+
+    return ProcessedLayer(output; labels, scale_mapping = merged_scale_mapping)
+end
diff --git a/src/transformations/highlight.jl b/src/transformations/highlight.jl
new file mode 100644
index 0000000000..cf4fde496b
--- /dev/null
+++ b/src/transformations/highlight.jl
@@ -0,0 +1,51 @@
+struct Highlight{F<:Function}
+    target::Vector{Union{Int,Symbol}}
+    predicate::F
+    repeat_facets::Bool
+end
+
+function Highlight(
+    target::Vector,
+    predicate;
+    repeat_facets = false
+)
+    return Highlight(
+        convert(Vector{Union{Int,Symbol}}, target),
+        predicate,
+        repeat_facets
+    )
+end
+
+function (highlight::Highlight)(p::ProcessedLayer)
+    primary = AlgebraOfGraphics.dictionary([(key == :color ? :group : key, value) for (key, value) in pairs(p.primary)
+        # should the data be repeated across all facets or not? maybe needs to be an option
+        if !(highlight.repeat_facets && key in (:layout, :row, :col))
+    ])
+    
+    grayed_out = ProcessedLayer(p; primary, attributes = merge(p.attributes, AlgebraOfGraphics.dictionary([:color => :gray80])))
+
+    function apply_predicate(target::Vector, predicate::Function, positional, named, scalar_primaries)
+        b = predicate([resolve_target(t, positional, named, scalar_primaries) for t in target]...)
+        b isa Bool || error("Highlighting predicate returned non-boolean value $b")
+        return b
+    end
+    
+    resolve_target(target::Int, positional, named, scalar_primaries) = positional[target]
+    resolve_target(target::Symbol, positional, named, scalar_primaries) = haskey(scalar_primaries, target) ? scalar_primaries[target] : named[target]
+
+
+    i::Int = 0
+    colored = AlgebraOfGraphics.filtermap(p) do positional, named
+        i += 1
+        scalar_primaries = map(val -> val[i], p.primary)
+        if apply_predicate(highlight.target, highlight.predicate, positional, named, scalar_primaries)
+            return positional, named
+        else
+            return nothing
+        end
+    end
+
+    return ProcessedLayers([grayed_out, colored])
+end
+
+highlight(pair::Pair; kwargs...) = AlgebraOfGraphics.transformation(Highlight(vcat(pair[1]), pair[2]; kwargs...))
\ No newline at end of file
diff --git a/test/analyses.jl b/test/analyses.jl
index ffbb7f1f65..4ee252adc1 100644
--- a/test/analyses.jl
+++ b/test/analyses.jl
@@ -660,3 +660,18 @@ end
     @test x̂[2] ≈ x̂2
     @test ŷ[2] ≈ ŷ2
 end
+
+@testset "aggregate fails with higher-dimensional result" begin
+    df = (x = [1, 1, 2, 2], y = [1, 2, 3, 4])
+
+    matrix_func = v -> reshape(v, 1, :)
+
+    layer = data(df) * mapping(:x, :y) * aggregate(2 => matrix_func)
+
+    @test_throws "Aggregation of positional argument 2 returned 2-dimensional arrays with size (1, 2). Only scalars or 1-dimensional vectors are supported." AlgebraOfGraphics.ProcessedLayer(layer)
+end
+
+@testset "aggregate fails for multiply assigned output slots" begin
+    @test_throws "Output slot 2 of `aggregate` was assigned multiple times." data((; x = [1, 1, 1], y = [1, 2, 3])) *
+        mapping(:x, :y) * aggregate(2 => mean, 2 => std) * visual(Errorbars) |> draw
+end
diff --git a/test/reference_tests.jl b/test/reference_tests.jl
index 662de15274..a84522b11f 100644
--- a/test/reference_tests.jl
+++ b/test/reference_tests.jl
@@ -1622,3 +1622,300 @@ reftest("all-missing groups") do
     draw!(f[1, 2], spec2)
     f
 end
+
+reftest("aggregate mean over x values") do
+    # Three groups with different numbers of points (4, 5, 6) and mean line going up and down
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            2.0, 2.2, 1.8, 2.0,  # mean = 2.0
+            4.8, 5.0, 5.2, 4.9, 5.1,  # mean = 5.0
+            2.7, 3.0, 3.3, 2.8, 3.2, 3.0,  # mean = 3.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y) * visual(Scatter, color = :gray)
+    layer_mean = data(df) * mapping(:x, :y) * aggregate(2 => mean) * visual(Lines, color = :red, linewidth = 3)
+    draw(layer_raw + layer_mean)
+end
+
+reftest("aggregate mean with layout faceting") do
+    # Two groups (A and B) with different linear relationships
+    df = (;
+        x = [
+            1, 1, 1, 2, 2, 2, 3, 3, 3,  # Group A x values
+            1, 1, 1, 2, 2, 2, 3, 3, 3,   # Group B x values
+        ],
+        y = [
+            1.8, 2.0, 2.2, 3.8, 4.0, 4.2, 5.7, 6.0, 6.3,  # Group A: y ≈ 2x (means: 2.0, 4.0, 6.0)
+            2.7, 3.0, 3.3, 5.7, 6.0, 6.3, 8.7, 9.0, 9.3,   # Group B: y ≈ 3x (means: 3.0, 6.0, 9.0)
+        ],
+        group = [
+            "A", "A", "A", "A", "A", "A", "A", "A", "A",
+            "B", "B", "B", "B", "B", "B", "B", "B", "B",
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y, layout = :group) * visual(Scatter, color = :gray)
+    layer_mean = data(df) * mapping(:x, :y, layout = :group) * aggregate(2 => mean) * visual(Lines, color = :red, linewidth = 3)
+    draw(layer_raw + layer_mean)
+end
+
+reftest("aggregate mean of x over y values") do
+    # Aggregate x values for each y group (horizontal aggregation) with zig-zag pattern
+    df = (;
+        y = [
+            1, 1, 1, 1,  # 4 points at y=1
+            2, 2, 2, 2, 2,  # 5 points at y=2
+            3, 3, 3, 3, 3, 3,  # 6 points at y=3
+        ],
+        x = [
+            1.8, 2.0, 2.2, 2.0,  # mean = 2.0
+            4.8, 5.0, 5.2, 4.9, 5.1,  # mean = 5.0
+            2.7, 3.0, 3.3, 2.8, 3.2, 3.0,  # mean = 3.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y) * visual(Scatter, color = :gray)
+    layer_mean = data(df) * mapping(:x, :y) * aggregate(1 => mean) * visual(Lines, color = :red, linewidth = 3)
+    draw(layer_raw + layer_mean)
+end
+
+reftest("aggregate mean with color aggregation") do
+    # Groups with different sizes (zig-zag pattern), color shows group size via length aggregation
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            1.8, 2.0, 2.2, 2.0,  # mean = 2.0
+            4.8, 5.0, 5.2, 4.9, 5.1,  # mean = 5.0
+            2.7, 3.0, 3.3, 2.8, 3.2, 3.0,  # mean = 3.0
+        ],
+        color = [
+            0.8, 1.0, 1.2, 1.0,  # mean = 1.0
+            1.8, 2.0, 2.2, 1.9, 2.1,  # mean = 2.0
+            2.7, 3.0, 3.3, 2.8, 3.2, 3.0,  # mean = 3.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y) * visual(Scatter, color = :gray)
+    layer_agg = data(df) * mapping(:x, :y, color = :color) *
+        aggregate(2 => mean, :color => length) *
+        visual(Scatter, markersize = 20, marker = :diamond, colormap = :viridis)
+    draw(layer_raw + layer_agg)
+end
+
+reftest("aggregate mean with missing values") do
+    # One group has a missing value - mean should return missing for that group
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2, one will be missing
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            1.8, 2.0, 2.2, 2.0,  # mean = 2.0
+            4.8, missing, 5.2, 4.9, 5.1,  # mean = missing (because one value is missing)
+            2.7, 3.0, 3.3, 2.8, 3.2, 3.0,  # mean = 3.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y) * visual(Scatter, color = :gray)
+    layer_mean = data(df) * mapping(:x, :y) * aggregate(2 => mean) * visual(Scatter, color = :blue, markersize = 20)
+    draw(layer_raw + layer_mean)
+end
+
+reftest("aggregate sum heatmap 2d") do
+    # 2x3 heatmap with one combination missing (x=2, y=2) to show gap
+    df = (;
+        x = [
+            1, 1, 1,  # (1,1): sum = 6.0
+            1, 1,  # (1,2): sum = 5.0
+            1, 1, 1,  # (1,3): sum = 9.0
+            2, 2,  # (2,1): sum = 7.0
+            # (2,2): missing - no data points
+            2, 2, 2, 2,  # (2,3): sum = 12.0
+        ],
+        y = [
+            1, 1, 1,
+            2, 2,
+            3, 3, 3,
+            1, 1,
+            3, 3, 3, 3,
+        ],
+        z = [
+            2.0, 2.0, 2.0,  # sum = 6.0
+            2.5, 2.5,  # sum = 5.0
+            3.0, 3.0, 3.0,  # sum = 9.0
+            3.0, 4.0,  # sum = 7.0
+            3.0, 3.0, 3.0, 3.0,  # sum = 12.0
+        ],
+    )
+    data(df) * mapping(:x, :y, :z) * aggregate(3 => sum) * visual(Heatmap) |> draw
+end
+
+reftest("aggregate extrema rangebars") do
+    # Extrema split into min and max for range bars with zig-zag pattern
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            1.5, 2.0, 2.5, 2.2,  # min = 1.5, max = 2.5
+            4.3, 5.0, 5.7, 4.8, 5.2,  # min = 4.3, max = 5.7
+            2.0, 3.0, 4.0, 2.5, 3.5, 3.2,  # min = 2.0, max = 4.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y) * visual(Scatter, color = :gray)
+    layer_range = data(df) * mapping(:x, :y) *
+        aggregate(2 => extrema => [first => 2, last => 3]) *
+        visual(Rangebars, color = :red, linewidth = 3)
+    draw(layer_raw + layer_range)
+end
+
+reftest("aggregate sum heatmap custom scale and label") do
+    df = (;
+        x = [
+            1, 1, 1,  # (1,1): sum = 6.0
+            1, 1,  # (1,2): sum = 5.0
+            1, 1, 1,  # (1,3): sum = 9.0
+            2, 2,  # (2,1): sum = 7.0
+            # (2,2): missing - no data points
+            2, 2, 2, 2,  # (2,3): sum = 12.0
+        ],
+        y = [
+            1, 1, 1,
+            2, 2,
+            3, 3, 3,
+            1, 1,
+            3, 3, 3, 3,
+        ],
+        z = [
+            2.0, 2.0, 2.0,  # sum = 6.0
+            2.5, 2.5,  # sum = 5.0
+            3.0, 3.0, 3.0,  # sum = 9.0
+            3.0, 4.0,  # sum = 7.0
+            3.0, 3.0, 3.0, 3.0,  # sum = 12.0
+        ],
+    )
+    layer = data(df) * mapping(:x, :y, :z) *
+        aggregate(3 => sum => rich("total ", rich("of z", font = :bold)) => scale(:color2)) *
+        visual(Heatmap)
+    draw(layer, scales(color2 = (; colormap = :Blues)))
+end
+
+reftest("aggregate extrema split custom labels and scale") do
+    # Split extrema: min becomes y position with label, max becomes color with custom scale
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            1.5, 2.0, 2.5, 2.2,  # min = 1.5, max = 2.5
+            4.3, 5.0, 5.7, 4.8, 5.2,  # min = 4.3, max = 5.7
+            2.0, 3.0, 4.0, 2.5, 3.5, 3.2,  # min = 2.0, max = 4.0
+        ],
+    )
+    layer_raw = data(df) * mapping(:x, :y => "") * visual(Scatter, color = :gray)
+    layer_agg = data(df) * mapping(:x, :y) *
+        aggregate(
+        2 => extrema => [
+            first => 2 => "Min",  # Lower bound as y coordinate with label "Min"
+            last => :color => "Max" => scale(:color2),  # Upper bound as color with label "Max" and custom scale
+        ]
+    ) *
+        visual(Scatter, markersize = 25)
+    draw(layer_raw + layer_agg, scales(color2 = (; colormap = :thermal)))
+end
+
+reftest("aggregate categorical from numerical") do
+    # Custom aggregation function that returns both category and mean value
+    function categorize_mean(values)
+        m = mean(values)
+        category = if m < 3.0
+            "low"
+        elseif m < 6.0
+            "mid"
+        else
+            "high"
+        end
+        return (category, m)
+    end
+
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1
+            2, 2, 2, 2, 2,  # 5 points at x=2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3
+        ],
+        y = [
+            1.8, 2.0, 2.2, 2.0,  # mean = 2.0 → ("low", 2.0)
+            4.8, 5.0, 5.2, 4.9, 5.1,  # mean = 5.0 → ("mid", 5.0)
+            6.7, 7.0, 7.3, 6.8, 7.2, 7.0,  # mean = 7.0 → ("high", 7.0)
+        ],
+    )
+    spec = data(df) * mapping(:x, :y) *
+        aggregate(2 => categorize_mean => [first => 2 => "category", verbatim ∘ last => :bar_labels]) *
+        visual(BarPlot, direction = :x)
+    draw(spec, scales(X = (; categories = ["low", "mid", "high"])))
+end
+
+reftest("aggregate vector valued") do
+    # Vector-valued aggregation: extract lower half of points in each group
+    function lower_half(values)
+        n = length(values)
+        n_lower = div(n, 2, RoundDown)
+        sorted = sort(values)
+        return sorted[1:n_lower]
+    end
+
+    df = (;
+        x = [
+            1, 1, 1, 1,  # 4 points at x=1 → lower 2
+            2, 2, 2, 2, 2,  # 5 points at x=2 → lower 2
+            3, 3, 3, 3, 3, 3,  # 6 points at x=3 → lower 3
+        ],
+        y = [
+            1.8, 2.0, 2.2, 2.4,  # lower half: [1.8, 2.0]
+            4.6, 4.8, 5.0, 5.2, 5.4,  # lower half: [4.6, 4.8]
+            2.5, 2.7, 3.0, 3.3, 3.5, 3.7,  # lower half: [2.5, 2.7, 3.0]
+        ],
+    )
+    # Show all points with larger markers
+    layer_all = data(df) * mapping(:x, :y) * visual(Scatter, markersize = 20, color = :gray80)
+    # Overlay lower half with smaller markers - should align perfectly with subset of gray points
+    layer_lower = data(df) * mapping(:x, :y) * aggregate(2 => lower_half) * visual(Scatter, markersize = 10, color = :red)
+    draw(layer_all + layer_lower)
+end
+
+reftest("aggregate no group") do
+    df = (; x = 1:5, y = 6:10)
+    spec = data(df) *
+        (
+        mapping(:x, :y) * visual(Scatter) +
+            mapping(:y) * aggregate(1 => collect ∘ extrema) * visual(HLines)
+    )
+    draw(spec)
+end
+
+reftest("aggregate tupleized multi input single output") do
+    df = (; a = [1, 1, 1, 2, 2], b = [0, -3, 4, 0, 1], c = [1, -4, 2, 0, 2])
+    data(df) * mapping(:a => nonnumeric, (:b, :c)) * aggregate(2 => (x -> mean([maximum(t) for t in x])) => "mean of maxs") * visual(Scatter) |> draw
+end
+
+reftest("two columns and assignment") do
+    df = (
+        x = [1, 1, 1, 1, 2, 2, 2, 2],
+        y = [1, 2, 3, 4, -2, 0, 4, 7],
+    )
+    data(df) *
+        mapping(:x => nonnumeric, :y) *
+        aggregate(2 => mean, 2 => std => 3) *
+        visual(Errorbars) |> draw
+end
diff --git a/test/reference_tests/aggregate categorical from numerical ref.png b/test/reference_tests/aggregate categorical from numerical ref.png
new file mode 100644
index 0000000000..5742b42a5c
Binary files /dev/null and b/test/reference_tests/aggregate categorical from numerical ref.png differ
diff --git a/test/reference_tests/aggregate extrema rangebars ref.png b/test/reference_tests/aggregate extrema rangebars ref.png
new file mode 100644
index 0000000000..616ec9a2bb
Binary files /dev/null and b/test/reference_tests/aggregate extrema rangebars ref.png differ
diff --git a/test/reference_tests/aggregate extrema split custom labels and scale ref.png b/test/reference_tests/aggregate extrema split custom labels and scale ref.png
new file mode 100644
index 0000000000..fe769a0450
Binary files /dev/null and b/test/reference_tests/aggregate extrema split custom labels and scale ref.png differ
diff --git a/test/reference_tests/aggregate mean of x over y values ref.png b/test/reference_tests/aggregate mean of x over y values ref.png
new file mode 100644
index 0000000000..ad010a06ed
Binary files /dev/null and b/test/reference_tests/aggregate mean of x over y values ref.png differ
diff --git a/test/reference_tests/aggregate mean over x values ref.png b/test/reference_tests/aggregate mean over x values ref.png
new file mode 100644
index 0000000000..c1be579c9c
Binary files /dev/null and b/test/reference_tests/aggregate mean over x values ref.png differ
diff --git a/test/reference_tests/aggregate mean with color aggregation ref.png b/test/reference_tests/aggregate mean with color aggregation ref.png
new file mode 100644
index 0000000000..514eac4ddd
Binary files /dev/null and b/test/reference_tests/aggregate mean with color aggregation ref.png differ
diff --git a/test/reference_tests/aggregate mean with layout faceting ref.png b/test/reference_tests/aggregate mean with layout faceting ref.png
new file mode 100644
index 0000000000..d6922a4599
Binary files /dev/null and b/test/reference_tests/aggregate mean with layout faceting ref.png differ
diff --git a/test/reference_tests/aggregate mean with missing values ref.png b/test/reference_tests/aggregate mean with missing values ref.png
new file mode 100644
index 0000000000..16abfd8f34
Binary files /dev/null and b/test/reference_tests/aggregate mean with missing values ref.png differ
diff --git a/test/reference_tests/aggregate no group ref.png b/test/reference_tests/aggregate no group ref.png
new file mode 100644
index 0000000000..b4b6c6ae18
Binary files /dev/null and b/test/reference_tests/aggregate no group ref.png differ
diff --git a/test/reference_tests/aggregate sum heatmap 2d ref.png b/test/reference_tests/aggregate sum heatmap 2d ref.png
new file mode 100644
index 0000000000..813c710410
Binary files /dev/null and b/test/reference_tests/aggregate sum heatmap 2d ref.png differ
diff --git a/test/reference_tests/aggregate sum heatmap custom scale and label ref.png b/test/reference_tests/aggregate sum heatmap custom scale and label ref.png
new file mode 100644
index 0000000000..7bd26f184e
Binary files /dev/null and b/test/reference_tests/aggregate sum heatmap custom scale and label ref.png differ
diff --git a/test/reference_tests/aggregate tupleized multi input single output ref.png b/test/reference_tests/aggregate tupleized multi input single output ref.png
new file mode 100644
index 0000000000..8528c7bd68
Binary files /dev/null and b/test/reference_tests/aggregate tupleized multi input single output ref.png differ
diff --git a/test/reference_tests/aggregate vector valued ref.png b/test/reference_tests/aggregate vector valued ref.png
new file mode 100644
index 0000000000..ba3ba4da6b
Binary files /dev/null and b/test/reference_tests/aggregate vector valued ref.png differ
diff --git a/test/reference_tests/two columns and assignment ref.png b/test/reference_tests/two columns and assignment ref.png
new file mode 100644
index 0000000000..3d84fa3694
Binary files /dev/null and b/test/reference_tests/two columns and assignment ref.png differ