diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml
index e3f0749..acce3e7 100644
--- a/.github/workflows/Tests.yml
+++ b/.github/workflows/Tests.yml
@@ -14,21 +14,7 @@ concurrency:
 
 jobs:
   tests:
-    name: "Tests"
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - "lts"
-          - "1"
-          - "pre"
-        os:
-          - ubuntu-latest
-          - macOS-latest
-          - windows-latest
-    uses: "SciML/.github/.github/workflows/tests.yml@v1"
+    uses: "SciML/.github/.github/workflows/grouped-tests.yml@v1"
     with:
-      julia-version: "${{ matrix.version }}"
-      os: "${{ matrix.os }}"
       coverage: true
     secrets: "inherit"
diff --git a/.gitignore b/.gitignore
index 73d4447..ea35671 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@
 *.jl.mem
 /Manifest.toml
 /docs/build/
-/docs/Manifest.toml
\ No newline at end of file
+/docs/Manifest.toml
+/test/qa/Manifest.toml
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index cb61929..2cd566d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,17 +7,20 @@ authors = ["ChrisRackauckas <contact@chrisrackauckas.com>"]
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 
 [compat]
-ExplicitImports = "1.14.0"
-JET = "0.9, 0.10, 0.11"
 PrecompileTools = "1.2"
+Random = "1.10"
+SafeTestsets = "0.1"
+SciMLTesting = "1.6"
+Statistics = "1.10"
+Test = "1.10"
 julia = "1.10"
 
 [extras]
-ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ExplicitImports", "JET", "Random", "Statistics", "Test"]
+test = ["Random", "SafeTestsets", "SciMLTesting", "Statistics", "Test"]
diff --git a/test/core_tests.jl b/test/core_tests.jl
new file mode 100644
index 0000000..456e725
--- /dev/null
+++ b/test/core_tests.jl
@@ -0,0 +1,135 @@
+using LightweightStats
+using Test
+
+using LightweightStats: mean, median, std, var, cov, cor, quantile, middle
+
+@testset "mean" begin
+    @test mean([1, 2, 3, 4, 5]) ≈ 3.0
+    @test mean([1.5, 2.5, 3.5]) ≈ 2.5
+    @test mean(Float32[1, 2, 3]) ≈ 2.0f0
+
+    @test mean(x -> x^2, [1, 2, 3]) ≈ 14 / 3
+
+    A = [1 2 3; 4 5 6]
+    @test mean(A) ≈ 3.5
+    @test mean(A; dims = 1) ≈ [2.5 3.5 4.5]
+    @test mean(A; dims = 2) ≈ reshape([2.0, 5.0], 2, 1)
+
+    @test_throws ArgumentError mean([])
+end
+
+@testset "median" begin
+    @test median([1, 2, 3, 4, 5]) == 3
+    @test median([1, 2, 3, 4]) == 2.5
+    @test median([3, 1, 2]) == 2
+    @test median([1]) == 1
+
+    A = [1 2 3; 4 5 6]
+    @test median(A) == 3.5
+    @test median(A; dims = 1) == [2.5 3.5 4.5]
+    @test median(A; dims = 2) == reshape([2.0, 5.0], 2, 1)
+
+    @test_throws ArgumentError median([])
+end
+
+@testset "var and std" begin
+    x = [1, 2, 3, 4, 5]
+    @test var(x) ≈ 2.5
+    @test var(x; corrected = false) ≈ 2.0
+    @test std(x) ≈ sqrt(2.5)
+    @test std(x; corrected = false) ≈ sqrt(2.0)
+
+    @test var(x; mean = 3) ≈ 2.5
+    @test std(x; mean = 3) ≈ sqrt(2.5)
+
+    A = [1 2 3; 4 5 6]
+    @test var(A) ≈ 3.5
+    @test std(A) ≈ sqrt(3.5)
+
+    @test isnan(var(Float64[]))
+end
+
+@testset "cov" begin
+    x = [1.0, 2.0, 3.0, 4.0, 5.0]
+    y = [2.0, 4.0, 6.0, 8.0, 10.0]
+
+    @test cov(x, y) ≈ 5.0
+    @test cov(x, y; corrected = false) ≈ 4.0
+    @test cov(x) ≈ var(x)
+
+    X = [1 2; 3 4; 5 6]
+    C = cov(X; dims = 1)
+    @test size(C) == (2, 2)
+    @test C[1, 1] ≈ var(X[:, 1])
+    @test C[2, 2] ≈ var(X[:, 2])
+    @test C[1, 2] ≈ C[2, 1]
+
+    C2 = cov(X; dims = 2)
+    @test size(C2) == (3, 3)
+
+    @test_throws DimensionMismatch cov([1, 2], [1, 2, 3])
+end
+
+@testset "cor" begin
+    x = [1.0, 2.0, 3.0, 4.0, 5.0]
+    y = [2.0, 4.0, 6.0, 8.0, 10.0]
+
+    @test cor(x, y) ≈ 1.0
+    @test cor(x, -y) ≈ -1.0
+
+    X = [1 2; 3 4; 5 6]
+    R = cor(X; dims = 1)
+    @test size(R) == (2, 2)
+    @test R[1, 1] ≈ 1.0
+    @test R[2, 2] ≈ 1.0
+    @test R[1, 2] ≈ R[2, 1]
+
+    @test isnan(cor([1, 1, 1], [1, 2, 3]))
+
+    @test_throws DimensionMismatch cor([1, 2], [1, 2, 3])
+end
+
+@testset "quantile" begin
+    v = [1, 2, 3, 4, 5]
+
+    @test quantile(v, 0.0) == 1
+    @test quantile(v, 0.25) ≈ 2.0
+    @test quantile(v, 0.5) ≈ 3.0
+    @test quantile(v, 0.75) ≈ 4.0
+    @test quantile(v, 1.0) == 5
+
+    @test quantile(v, [0.25, 0.5, 0.75]) ≈ [2.0, 3.0, 4.0]
+
+    @test_throws ArgumentError quantile(v, -0.1)
+    @test_throws ArgumentError quantile(v, 1.1)
+    @test_throws ArgumentError quantile([], 0.5)
+end
+
+@testset "middle" begin
+    @test middle(1, 5) == 3
+    @test middle(1.0, 5.0) == 3.0
+    @test middle([1, 2, 3, 4, 5]) == 3
+    @test middle([5, 1, 3]) == 3
+    @test middle(10) == 10
+
+    @test_throws ArgumentError middle([])
+end
+
+@testset "Type stability" begin
+    @test typeof(mean(Int[1, 2, 3])) == Float64
+    @test typeof(mean(Float32[1, 2, 3])) == Float32
+    @test typeof(std(Int[1, 2, 3])) == Float64
+    @test typeof(median(Int[1, 2, 3])) == Int
+    @test typeof(median(Float64[1, 2, 3])) == Float64
+end
+
+@testset "Edge cases" begin
+    @test mean([42]) == 42
+    @test median([42]) == 42
+    @test var([42]) |> isnan
+    @test std([42]) |> isnan
+
+    @test mean([1, 1, 1, 1]) == 1
+    @test var([1, 1, 1, 1]; corrected = false) == 0
+    @test std([1, 1, 1, 1]; corrected = false) == 0
+end
diff --git a/test/explicit_imports_tests.jl b/test/explicit_imports_tests.jl
deleted file mode 100644
index dfe8a51..0000000
--- a/test/explicit_imports_tests.jl
+++ /dev/null
@@ -1,8 +0,0 @@
-using ExplicitImports
-using LightweightStats
-using Test
-
-@testset "ExplicitImports" begin
-    @test check_no_implicit_imports(LightweightStats) === nothing
-    @test check_no_stale_explicit_imports(LightweightStats) === nothing
-end
diff --git a/test/jet_tests.jl b/test/jet_tests.jl
deleted file mode 100644
index f9f0af0..0000000
--- a/test/jet_tests.jl
+++ /dev/null
@@ -1,126 +0,0 @@
-using LightweightStats
-using Test
-
-# JET tests are optional and only run on stable Julia releases 1.11+
-# JET has strict Julia version requirements and tight coupling with the compiler:
-# - JET 0.9.x works with Julia 1.10/1.11
-# - JET 0.10.x/0.11.x works with Julia 1.12+
-# We skip JET tests on:
-# - Julia < 1.11 (LTS may have JET API differences)
-# - Julia pre-release/nightly (may not have JET support yet)
-const IS_STABLE_JULIA = VERSION >= v"1.11" && isempty(VERSION.prerelease)
-const JET_AVAILABLE = IS_STABLE_JULIA && try
-    @eval using JET
-    true
-catch e
-    @info "JET not available: $e"
-    false
-end
-
-@testset "JET static analysis" begin
-    if !JET_AVAILABLE
-        @info "JET tests skipped on Julia $(VERSION)"
-        @test_skip true  # Mark test as skipped
-        return
-    end
-    @testset "JET error analysis" begin
-        # Test key entry points for static errors using report_call
-        rep = JET.report_call(mean, (Vector{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(mean, (Vector{Int},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(median, (Vector{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(var, (Vector{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(std, (Vector{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cov, (Vector{Float64}, Vector{Float64}))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cov, (Matrix{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cor, (Vector{Float64}, Vector{Float64}))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cor, (Matrix{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(quantile, (Vector{Float64}, Float64))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(middle, (Vector{Float64},))
-        @test length(JET.get_reports(rep)) == 0
-    end
-
-    @testset "JET optimization analysis" begin
-        # Test key entry points for type stability using report_opt
-        # Use target_modules to filter to LightweightStats only
-        # (Base.mapslices has known runtime dispatches that are not our concern)
-        rep = JET.report_opt(mean, (Vector{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(median, (Vector{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(var, (Vector{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(std, (Vector{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(cov, (Vector{Float64}, Vector{Float64}); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(cov, (Matrix{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(cor, (Vector{Float64}, Vector{Float64}); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(cor, (Matrix{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(quantile, (Vector{Float64}, Float64); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_opt(middle, (Vector{Float64},); target_modules = (LightweightStats,))
-        @test length(JET.get_reports(rep)) == 0
-    end
-
-    @testset "JET analysis with different numeric types" begin
-        # Float32 type preservation
-        rep = JET.report_call(mean, (Vector{Float32},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(var, (Vector{Float32},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cov, (Matrix{Float32},))
-        @test length(JET.get_reports(rep)) == 0
-
-        # BigFloat support
-        rep = JET.report_call(mean, (Vector{BigFloat},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(var, (Vector{BigFloat},))
-        @test length(JET.get_reports(rep)) == 0
-    end
-
-    @testset "JET analysis with complex numbers" begin
-        rep = JET.report_call(mean, (Vector{ComplexF64},))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cov, (Vector{ComplexF64}, Vector{ComplexF64}))
-        @test length(JET.get_reports(rep)) == 0
-
-        rep = JET.report_call(cov, (Matrix{ComplexF64},))
-        @test length(JET.get_reports(rep)) == 0
-    end
-end
diff --git a/test/qa/Project.toml b/test/qa/Project.toml
new file mode 100644
index 0000000..686d1af
--- /dev/null
+++ b/test/qa/Project.toml
@@ -0,0 +1,16 @@
+[deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+LightweightStats = "6e4b0b68-1f2b-4a7e-8f0e-d3e9f3d5c8a1"
+SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[sources]
+LightweightStats = {path = "../.."}
+
+[compat]
+Aqua = "0.8"
+JET = "0.9, 0.10, 0.11"
+SciMLTesting = "1.6"
+Test = "1"
+julia = "1.10"
diff --git a/test/qa/qa.jl b/test/qa/qa.jl
new file mode 100644
index 0000000..951aa91
--- /dev/null
+++ b/test/qa/qa.jl
@@ -0,0 +1,3 @@
+using SciMLTesting, LightweightStats, JET, Test
+
+run_qa(LightweightStats; explicit_imports = true)
diff --git a/test/runtests.jl b/test/runtests.jl
index fda0661..a18a7cc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,147 +1,2 @@
-using LightweightStats
-using Test
-using Random
-
-# Import functions from LightweightStats for cleaner test code
-using LightweightStats: mean, median, std, var, cov, cor, quantile, middle
-
-@testset "LightweightStats.jl" begin
-    # Include regression tests against Statistics.jl
-    include("regression_tests.jl")
-    # Include interface compatibility tests
-    include("interface_tests.jl")
-    # Include JET static analysis tests
-    include("jet_tests.jl")
-    # Include explicit imports tests
-    include("explicit_imports_tests.jl")
-    @testset "mean" begin
-        @test mean([1, 2, 3, 4, 5]) ≈ 3.0
-        @test mean([1.5, 2.5, 3.5]) ≈ 2.5
-        @test mean(Float32[1, 2, 3]) ≈ 2.0f0
-
-        @test mean(x -> x^2, [1, 2, 3]) ≈ 14 / 3
-
-        A = [1 2 3; 4 5 6]
-        @test mean(A) ≈ 3.5
-        @test mean(A; dims = 1) ≈ [2.5 3.5 4.5]
-        @test mean(A; dims = 2) ≈ reshape([2.0, 5.0], 2, 1)
-
-        @test_throws ArgumentError mean([])
-    end
-
-    @testset "median" begin
-        @test median([1, 2, 3, 4, 5]) == 3
-        @test median([1, 2, 3, 4]) == 2.5
-        @test median([3, 1, 2]) == 2
-        @test median([1]) == 1
-
-        A = [1 2 3; 4 5 6]
-        @test median(A) == 3.5
-        @test median(A; dims = 1) == [2.5 3.5 4.5]
-        @test median(A; dims = 2) == reshape([2.0, 5.0], 2, 1)
-
-        @test_throws ArgumentError median([])
-    end
-
-    @testset "var and std" begin
-        x = [1, 2, 3, 4, 5]
-        @test var(x) ≈ 2.5
-        @test var(x; corrected = false) ≈ 2.0
-        @test std(x) ≈ sqrt(2.5)
-        @test std(x; corrected = false) ≈ sqrt(2.0)
-
-        @test var(x; mean = 3) ≈ 2.5
-        @test std(x; mean = 3) ≈ sqrt(2.5)
-
-        A = [1 2 3; 4 5 6]
-        @test var(A) ≈ 3.5
-        @test std(A) ≈ sqrt(3.5)
-
-        @test isnan(var(Float64[]))
-    end
-
-    @testset "cov" begin
-        x = [1.0, 2.0, 3.0, 4.0, 5.0]
-        y = [2.0, 4.0, 6.0, 8.0, 10.0]
-
-        @test cov(x, y) ≈ 5.0
-        @test cov(x, y; corrected = false) ≈ 4.0
-        @test cov(x) ≈ var(x)
-
-        X = [1 2; 3 4; 5 6]
-        C = cov(X; dims = 1)
-        @test size(C) == (2, 2)
-        @test C[1, 1] ≈ var(X[:, 1])
-        @test C[2, 2] ≈ var(X[:, 2])
-        @test C[1, 2] ≈ C[2, 1]
-
-        C2 = cov(X; dims = 2)
-        @test size(C2) == (3, 3)
-
-        @test_throws DimensionMismatch cov([1, 2], [1, 2, 3])
-    end
-
-    @testset "cor" begin
-        x = [1.0, 2.0, 3.0, 4.0, 5.0]
-        y = [2.0, 4.0, 6.0, 8.0, 10.0]
-
-        @test cor(x, y) ≈ 1.0
-        @test cor(x, -y) ≈ -1.0
-
-        X = [1 2; 3 4; 5 6]
-        R = cor(X; dims = 1)
-        @test size(R) == (2, 2)
-        @test R[1, 1] ≈ 1.0
-        @test R[2, 2] ≈ 1.0
-        @test R[1, 2] ≈ R[2, 1]
-
-        @test isnan(cor([1, 1, 1], [1, 2, 3]))
-
-        @test_throws DimensionMismatch cor([1, 2], [1, 2, 3])
-    end
-
-    @testset "quantile" begin
-        v = [1, 2, 3, 4, 5]
-
-        @test quantile(v, 0.0) == 1
-        @test quantile(v, 0.25) ≈ 2.0
-        @test quantile(v, 0.5) ≈ 3.0
-        @test quantile(v, 0.75) ≈ 4.0
-        @test quantile(v, 1.0) == 5
-
-        @test quantile(v, [0.25, 0.5, 0.75]) ≈ [2.0, 3.0, 4.0]
-
-        @test_throws ArgumentError quantile(v, -0.1)
-        @test_throws ArgumentError quantile(v, 1.1)
-        @test_throws ArgumentError quantile([], 0.5)
-    end
-
-    @testset "middle" begin
-        @test middle(1, 5) == 3
-        @test middle(1.0, 5.0) == 3.0
-        @test middle([1, 2, 3, 4, 5]) == 3
-        @test middle([5, 1, 3]) == 3
-        @test middle(10) == 10
-
-        @test_throws ArgumentError middle([])
-    end
-
-    @testset "Type stability" begin
-        @test typeof(mean(Int[1, 2, 3])) == Float64
-        @test typeof(mean(Float32[1, 2, 3])) == Float32
-        @test typeof(std(Int[1, 2, 3])) == Float64
-        @test typeof(median(Int[1, 2, 3])) == Int
-        @test typeof(median(Float64[1, 2, 3])) == Float64
-    end
-
-    @testset "Edge cases" begin
-        @test mean([42]) == 42
-        @test median([42]) == 42
-        @test var([42]) |> isnan
-        @test std([42]) |> isnan
-
-        @test mean([1, 1, 1, 1]) == 1
-        @test var([1, 1, 1, 1]; corrected = false) == 0
-        @test std([1, 1, 1, 1]; corrected = false) == 0
-    end
-end
+using SciMLTesting
+run_tests()
diff --git a/test/test_groups.toml b/test/test_groups.toml
new file mode 100644
index 0000000..f783870
--- /dev/null
+++ b/test/test_groups.toml
@@ -0,0 +1,18 @@
+# Root test-group matrix for LightweightStats.jl, consumed by the reusable
+# grouped-tests.yml workflow and by SciMLTesting's run_tests() folder-discovery
+# mode in runtests.jl. Each table key is a GROUP value passed to runtests.jl.
+#
+#   "Core" -> the top-level test/*.jl files (the functional + regression tests).
+#             Runs on lts/1/pre across the three OSes.
+#   "QA"   -> the test/qa/ folder (SciMLTesting run_qa: Aqua + ExplicitImports +
+#             JET), which is brittle on prereleases, so it runs only on lts/1.
+#             QA is always excluded from a local GROUP=All aggregate run.
+
+[Core]
+versions = ["lts", "1", "pre"]
+os = ["ubuntu-latest", "macOS-latest", "windows-latest"]
+
+[QA]
+in_all = false
+versions = ["lts", "1"]
+os = ["ubuntu-latest"]