diff --git a/Project.toml b/Project.toml
index cdecf98fb..7c63d7a84 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 
 [compat]
+Aqua = "0.6.5"
 DataAPI = "1"
 DataStructures = "0.10, 0.11, 0.12, 0.13, 0.14, 0.17, 0.18"
 LogExpFunctions = "0.3"
@@ -26,6 +27,7 @@ StatsAPI = "1.2"
 julia = "1"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
@@ -34,4 +36,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "Dates", "DelimitedFiles", "OffsetArrays", "StableRNGs", "Test"]
+test = ["Aqua", "BenchmarkTools", "Dates", "DelimitedFiles", "OffsetArrays", "StableRNGs", "Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 525749756..7fe605752 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,7 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Documenter = "0.27"
diff --git a/src/statmodels.jl b/src/statmodels.jl
index 291b6e2e7..89ce7b99a 100644
--- a/src/statmodels.jl
+++ b/src/statmodels.jl
@@ -55,7 +55,7 @@ end
 Show a p-value using 6 characters, either using the standard 0.XXXX
 representation or as <Xe-YY.
 """
-struct PValue <: Real
+struct PValue
     v::Real
     function PValue(v::Real)
         0 <= v <= 1 || isnan(v) || error("p-values must be in [0; 1]")
@@ -76,36 +76,22 @@ function show(io::IO, pv::PValue)
 end
 
 """Show a test statistic using 2 decimal digits"""
-struct TestStat <: Real
+struct TestStat
     v::Real
 end
-
-show(io::IO, x::TestStat) = @printf(io, "%.2f", x.v)
 TestStat(x::TestStat) = x
 
-float(x::Union{TestStat, PValue}) = float(x.v)
+show(io::IO, x::TestStat) = @printf(io, "%.2f", x.v)
 
-for op in [:(==), :<, :≤, :(isless), :(isequal)] # isless and < to place nice with NaN
-    @eval begin
-        Base.$op(x::Union{TestStat, PValue}, y::Real) = $op(x.v, y)
-        Base.$op(y::Real, x::Union{TestStat, PValue}) = $op(y, x.v)
-        Base.$op(x1::Union{TestStat, PValue}, x2::Union{TestStat, PValue}) = $op(x1.v, x2.v)
-    end
+# Better alignment when printing tables
+# Adapted from alignment(::IO, ::Real) (https://github.com/JuliaLang/julia/blob/0da46e25c865a390b5c2de20c2d40afb41fcac0a/base/show.jl#L3022-L3027)
+function Base.alignment(io::IO, x::Union{PValue,TestStat})
+    s = sprint(show, x)
+    m = match(r"^(.*?)((?:[\.eEfF].*)?)$", s)
+    m === nothing ? (textwidth(s), 0) :
+                    (textwidth(m.captures[1]), textwidth(m.captures[2]))
 end
 
-Base.hash(x::Union{TestStat, PValue}, h::UInt) = hash(x.v, h)
-
-# necessary to avoid a method ambiguity with isless(::TestStat, NaN)
-Base.isless(x::Union{TestStat, PValue}, y::AbstractFloat) = isless(x.v, y)
-Base.isless(y::AbstractFloat, x::Union{TestStat, PValue},) = isless(y, x.v)
-Base.isequal(y::AbstractFloat, x::Union{TestStat, PValue}) = isequal(y, x.v)
-Base.isequal(x::Union{TestStat, PValue}, y::AbstractFloat) = isequal(x.v, y)
-
-Base.isapprox(x::Union{TestStat, PValue}, y::Real; kwargs...) = isapprox(x.v, y; kwargs...)
-Base.isapprox(y::Real, x::Union{TestStat, PValue}; kwargs...) = isapprox(y, x.v; kwargs...)
-Base.isapprox(x1::Union{TestStat, PValue}, x2::Union{TestStat, PValue}; kwargs...) = isapprox(x1.v, x2.v; kwargs...)
-
-
 """Wrap a string so that show omits quotes"""
 struct NoQuote
     s::String
diff --git a/test/ambiguous.jl b/test/ambiguous.jl
deleted file mode 100644
index 59014ca0f..000000000
--- a/test/ambiguous.jl
+++ /dev/null
@@ -1,8 +0,0 @@
-using StatsBase
-using Test
-using Statistics
-
-@testset "Ambiguities" begin
-    # Ambiguities with Base, Core, and stdlib Statistics introduced by this package
-    @test_broken isempty(detect_ambiguities(StatsBase, Base, Core, Statistics))
-end
diff --git a/test/aqua.jl b/test/aqua.jl
new file mode 100644
index 000000000..fbddd7d60
--- /dev/null
+++ b/test/aqua.jl
@@ -0,0 +1,16 @@
+using Aqua
+using StatsBase
+using Test
+
+@testset "Aqua" begin
+    Aqua.test_all(
+        StatsBase;
+        # Ignore an ambiguity introduced by SortingAlgorithms
+        # Disable test on Julia < 1.6 due to ambiguities in ChainRulesCore
+        # (fixed in more recent versions not supported by Julia < 1.6)
+        ambiguities=VERSION < v"1.6" ? false : (exclude = [sort!],),
+        # `describe`, `pairwise`, `pairwise!`, and `weights` are defined in StatsAPI but owned by StasBase
+        # `quantile(::AbstractVector{<:Real})` is a deprecated method that should be removed
+        piracy=(treat_as_own = [describe, pairwise, pairwise!, weights, quantile],),
+    ) 
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 94986d240..fe4828a8a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,8 +2,9 @@ using StatsBase
 using Dates
 using LinearAlgebra
 using Random
+using Statistics
 
-tests = ["ambiguous",
+tests = ["aqua",
          "weights",
          "moments",
          "scalarstats",
diff --git a/test/statmodels.jl b/test/statmodels.jl
index 778ea2e39..ec45f5b3c 100644
--- a/test/statmodels.jl
+++ b/test/statmodels.jl
@@ -98,29 +98,6 @@ end
 @test sprint(show, TestStat(1e-5)) == "0.00"
 @test sprint(show, TestStat(π)) == "3.14"
 
-@testset "Union{PValue, TestStat} is Real" begin
-    vals = [0.0, Rational(1,3), NaN]
-    for T in [PValue, TestStat],
-        f in (==, <, ≤, >, ≥, isless, isequal),
-        lhs in vals, rhs in vals
-        # make sure that T behaves like a Real,
-        # regardless of whether it's on the LHS, RHS or both
-        @test f(T(lhs), T(rhs)) == f(lhs, rhs)
-        @test f(lhs, T(rhs)) == f(lhs, rhs)
-        @test f(T(lhs), rhs) == f(lhs, rhs)
-    end
-
-    # the (approximate) equality operators get a bit more attention
-    for T in [PValue, TestStat]
-        @test T(Rational(1,3)) ≈ T(1/3)
-        @test Rational(1,3) ≈ T(1/3) atol=0.01
-        @test T(Rational(1,3)) isa Real
-        @test T(T(0.05)) === T(0.05)
-        @test hash(T(0.05)) == hash(0.05)
-        @test hash(T(0.05), UInt(42)) == hash(0.05, UInt(42))
-    end
-end
-
 @test sprint(showerror, ConvergenceException(10)) == "failure to converge after 10 iterations."
 
 @test sprint(showerror, ConvergenceException(10, 0.2, 0.1)) ==