diff --git a/Project.toml b/Project.toml
index 69b4e0d..f115137 100644
--- a/Project.toml
+++ b/Project.toml
@@ -28,6 +28,7 @@ HypergeometricFunctions = "0.3.10"
 InverseFunctions = "0.1"
 IrrationalConstants = "0.1, 0.2"
 JET = "0.9, 0.10, 0.11"
+TestItemRunner = "1"
 LogExpFunctions = "0.3.26"
 Random = "<0.0.1, 1"
 Reexport = "1"
@@ -46,6 +47,7 @@ InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Aqua", "ChainRulesCore", "ChainRulesTestUtils", "ExplicitImports", "ForwardDiff", "InverseFunctions", "JET", "Random", "Test"]
+test = ["Aqua", "ChainRulesCore", "ChainRulesTestUtils", "ExplicitImports", "ForwardDiff", "InverseFunctions", "JET", "Random", "Test", "TestItemRunner"]
diff --git a/test/chainrules.jl b/test/chainrules.jl
index d59b5ea..f9a6436 100644
--- a/test/chainrules.jl
+++ b/test/chainrules.jl
@@ -1,53 +1,55 @@
-using StatsFuns, Test
-using ChainRulesCore
-using ChainRulesTestUtils
-using Random
-
-@testset "chainrules" begin
-    x = exp(randn())
-    y = exp(randn())
-    z = logistic(randn())
-    test_frule(betalogpdf, x, y, z)
-    test_rrule(betalogpdf, x, y, z)
-
-    x = exp(randn())
-    y = exp(randn())
-    z = exp(randn())
-    test_frule(gammalogpdf, x, y, z)
-    test_rrule(gammalogpdf, x, y, z)
-
-    x = exp(randn())
-    y = exp(randn())
-    test_frule(chisqlogpdf, x, y)
-    test_rrule(chisqlogpdf, x, y)
-
-    x = exp(randn())
-    y = exp(randn())
-    z = exp(randn())
-    test_frule(fdistlogpdf, x, y, z)
-    test_rrule(fdistlogpdf, x, y, z)
-
-    x = exp(randn())
-    y = randn()
-    test_frule(tdistlogpdf, x, y)
-    test_rrule(tdistlogpdf, x, y)
-
-    x = rand(1:100)
-    y = logistic(randn())
-    z = rand(1:x)
-    test_frule(binomlogpdf, x, y, z)
-    test_rrule(binomlogpdf, x, y, z)
-
-    x = exp(randn())
-    y = rand(1:100)
-    test_frule(poislogpdf, x, y)
-    test_rrule(poislogpdf, x, y)
-
-    # test special case λ = 0
-    _, pb = rrule(poislogpdf, 0.0, 0)
-    _, x̄1, _ = pb(1)
-    @test x̄1 == -1
-    _, pb = rrule(poislogpdf, 0.0, 1)
-    _, x̄1, _ = pb(1)
-    @test x̄1 == Inf
+@testitem "ChainRules" begin
+    using StatsFuns, Test
+    using ChainRulesCore
+    using ChainRulesTestUtils
+    using Random
+
+    @testset "chainrules" begin
+        x = exp(randn())
+        y = exp(randn())
+        z = logistic(randn())
+        test_frule(betalogpdf, x, y, z)
+        test_rrule(betalogpdf, x, y, z)
+
+        x = exp(randn())
+        y = exp(randn())
+        z = exp(randn())
+        test_frule(gammalogpdf, x, y, z)
+        test_rrule(gammalogpdf, x, y, z)
+
+        x = exp(randn())
+        y = exp(randn())
+        test_frule(chisqlogpdf, x, y)
+        test_rrule(chisqlogpdf, x, y)
+
+        x = exp(randn())
+        y = exp(randn())
+        z = exp(randn())
+        test_frule(fdistlogpdf, x, y, z)
+        test_rrule(fdistlogpdf, x, y, z)
+
+        x = exp(randn())
+        y = randn()
+        test_frule(tdistlogpdf, x, y)
+        test_rrule(tdistlogpdf, x, y)
+
+        x = rand(1:100)
+        y = logistic(randn())
+        z = rand(1:x)
+        test_frule(binomlogpdf, x, y, z)
+        test_rrule(binomlogpdf, x, y, z)
+
+        x = exp(randn())
+        y = rand(1:100)
+        test_frule(poislogpdf, x, y)
+        test_rrule(poislogpdf, x, y)
+
+        # test special case λ = 0
+        _, pb = rrule(poislogpdf, 0.0, 0)
+        _, x̄1, _ = pb(1)
+        @test x̄1 == -1
+        _, pb = rrule(poislogpdf, 0.0, 1)
+        _, x̄1, _ = pb(1)
+        @test x̄1 == Inf
+    end
 end
diff --git a/test/generic.jl b/test/generic.jl
index afe45e9..0acc338 100644
--- a/test/generic.jl
+++ b/test/generic.jl
@@ -1,50 +1,50 @@
-using StatsFuns
-using ForwardDiff: Dual
-using Test
-
-include("utils.jl")
-
-function check_dualx(statsfun, params, a, isprob, rtol)
-    v = @inferred(statsfun(params..., a))
-    rv = @inferred(statsfun(params..., Dual(a))).value
-    @test v isa float(Base.promote_typeof(params..., a))
-    @test rv isa float(Base.promote_typeof(params..., a))
-    return if isprob
-        @test v ≈ rv rtol = rtol nans = true
-    else
-        @test v ≈ rv atol = rtol rtol = rtol nans = true
+@testitem "Generic" setup=[Utils] begin
+    using StatsFuns
+    using ForwardDiff: Dual
+    using Test
+
+    _default_rtol = Utils._default_rtol
+
+    function check_dualx(statsfun, params, a, isprob, rtol)
+        v = @inferred(statsfun(params..., a))
+        rv = @inferred(statsfun(params..., Dual(a))).value
+        @test v isa float(Base.promote_typeof(params..., a))
+        @test rv isa float(Base.promote_typeof(params..., a))
+        return if isprob
+            @test v ≈ rv rtol = rtol nans = true
+        else
+            @test v ≈ rv atol = rtol rtol = rtol nans = true
+        end
     end
-end
 
-function genericcomp(basename::String, params, X::AbstractArray, rtol = _default_rtol(params, X))
-    if isdefined(@__MODULE__, Symbol(basename, :pdf))
-        stats_pdf = getproperty(@__MODULE__, Symbol(basename, :pdf))
-        @testset "pdf with params=$params, x=$x" for x in X
-            check_dualx(stats_pdf, params, x, true, rtol)
+    function genericcomp(basename::String, params, X::AbstractArray, rtol = _default_rtol(params, X))
+        if isdefined(@__MODULE__, Symbol(basename, :pdf))
+            stats_pdf = getproperty(@__MODULE__, Symbol(basename, :pdf))
+            @testset "pdf with params=$params, x=$x" for x in X
+                check_dualx(stats_pdf, params, x, true, rtol)
+            end
         end
-    end
 
-    if isdefined(@__MODULE__, Symbol(basename, :logpdf))
-        stats_logpdf = getproperty(@__MODULE__, Symbol(basename, :logpdf))
-        @testset "logpdf with params=$params, x=$x" for x in X
-            check_dualx(stats_logpdf, params, x, false, rtol)
+        if isdefined(@__MODULE__, Symbol(basename, :logpdf))
+            stats_logpdf = getproperty(@__MODULE__, Symbol(basename, :logpdf))
+            @testset "logpdf with params=$params, x=$x" for x in X
+                check_dualx(stats_logpdf, params, x, false, rtol)
+            end
         end
-    end
 
-    return nothing
-end
+        return nothing
+    end
 
-function genericcomp_tests(basename::String, configs)
-    println("\ttesting $basename ...")
-    for (params, data) in configs
-        genericcomp(basename, params, data)
+    function genericcomp_tests(basename::String, configs)
+        println("\ttesting $basename ...")
+        for (params, data) in configs
+            genericcomp(basename, params, data)
+        end
+        return
     end
-    return
-end
 
-### Test cases
+    ### Test cases
 
-@testset "Generic" begin
     genericcomp_tests(
         "beta", [
             ((1.0, 1.0), 0.01:0.01:0.99),
diff --git a/test/inverse.jl b/test/inverse.jl
index 1c7389a..b5e641f 100644
--- a/test/inverse.jl
+++ b/test/inverse.jl
@@ -1,11 +1,13 @@
-using StatsFuns
-using Test
-using InverseFunctions
+@testitem "Inverse" begin
+    using StatsFuns
+    using Test
+    using InverseFunctions
 
-@testset "inverse" begin
-    x = 0.7
-    InverseFunctions.test_inverse(norminvcdf, x)
-    InverseFunctions.test_inverse(norminvccdf, x)
-    InverseFunctions.test_inverse(norminvlogcdf, log(x))
-    InverseFunctions.test_inverse(norminvlogccdf, log(x))
+    @testset "inverse" begin
+        x = 0.7
+        InverseFunctions.test_inverse(norminvcdf, x)
+        InverseFunctions.test_inverse(norminvccdf, x)
+        InverseFunctions.test_inverse(norminvlogcdf, log(x))
+        InverseFunctions.test_inverse(norminvlogccdf, log(x))
+    end
 end
diff --git a/test/misc.jl b/test/misc.jl
index 40fa1c0..378b454 100644
--- a/test/misc.jl
+++ b/test/misc.jl
@@ -1,98 +1,100 @@
-using SpecialFunctions, StatsFuns
-using Test
+@testitem "Misc" begin
+    using SpecialFunctions, StatsFuns
+    using Test
 
-@testset "logmvgamma" begin
-    @testset "type behavior" for eltya in (Float32, Float64)
-        p = rand(1:50)
-        a = rand(eltya)
-        # add p since loggamma is only define for positive arguments
-        @test typeof(logmvgamma(p, a + p)) == eltya
-    end
-
-    @testset "consistent with loggamma" for eltya in (Float32, Float64)
-        #  Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc
-        a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments
-        @test logmvgamma(1, a) ≈ loggamma(a)
-        @test logmvgamma(2, a) ≈ eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5))
-        @test logmvgamma(3, a) ≈ eltya((3 / 2) * logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
-    end
+    @testset "logmvgamma" begin
+        @testset "type behavior" for eltya in (Float32, Float64)
+            p = rand(1:50)
+            a = rand(eltya)
+            # add p since loggamma is only define for positive arguments
+            @test typeof(logmvgamma(p, a + p)) == eltya
+        end
 
-    @testset "consistent with itself" for eltya in (Float32, Float64)
-        #  Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5)
-        for p in 1:50
-            a = rand(eltya) + p # add p since loggamma is only define for positive arguments
-            @test logmvgamma(p, a) ≈ eltya((p / 2 - 1 / 2) * logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
+        @testset "consistent with loggamma" for eltya in (Float32, Float64)
+            #  Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc
+            a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments
+            @test logmvgamma(1, a) ≈ loggamma(a)
+            @test logmvgamma(2, a) ≈ eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5))
+            @test logmvgamma(3, a) ≈ eltya((3 / 2) * logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a))
         end
-    end
-end
 
-@testset "logmvbeta" begin
-    @testset "symmetry" for eltya in (Float32, Float64)
-        for eltyb in (Float32, Float64)
-            #  Bᵢ(a, b) = Bᵢ(b, a)
+        @testset "consistent with itself" for eltya in (Float32, Float64)
+            #  Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5)
             for p in 1:50
-                a = rand(eltya) + p
-                b = rand(eltyb) + p
-                @test logmvbeta(p, a, b) ≈ logmvbeta(p, b, a)
+                a = rand(eltya) + p # add p since loggamma is only define for positive arguments
+                @test logmvgamma(p, a) ≈ eltya((p / 2 - 1 / 2) * logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5))
             end
         end
     end
 
-    @testset "consistent with logbeta" for eltya in (Float32, Float64)
-        for eltyb in (Float32, Float64)
-            #  B₁(a, b) = B(a, b)
-            a = rand(eltya)
-            b = rand(eltyb)
-            @test logmvbeta(1, a, b) ≈ logbeta(a, b)
+    @testset "logmvbeta" begin
+        @testset "symmetry" for eltya in (Float32, Float64)
+            for eltyb in (Float32, Float64)
+                #  Bᵢ(a, b) = Bᵢ(b, a)
+                for p in 1:50
+                    a = rand(eltya) + p
+                    b = rand(eltyb) + p
+                    @test logmvbeta(p, a, b) ≈ logmvbeta(p, b, a)
+                end
+            end
         end
-    end
 
-    @testset "type promotion behaves" for eltya in (Float32, Float64)
-        for eltyb in (Float32, Float64)
-            a = rand(eltya)
-            b = rand(eltyb)
-            T = Base.promote_type(eltya, eltyb)
-            @test typeof(logmvbeta(1, a, b)) == T
+        @testset "consistent with logbeta" for eltya in (Float32, Float64)
+            for eltyb in (Float32, Float64)
+                #  B₁(a, b) = B(a, b)
+                a = rand(eltya)
+                b = rand(eltyb)
+                @test logmvbeta(1, a, b) ≈ logbeta(a, b)
+            end
+        end
+
+        @testset "type promotion behaves" for eltya in (Float32, Float64)
+            for eltyb in (Float32, Float64)
+                a = rand(eltya)
+                b = rand(eltyb)
+                T = Base.promote_type(eltya, eltyb)
+                @test typeof(logmvbeta(1, a, b)) == T
+            end
         end
     end
-end
 
-# https://github.com/JuliaStats/StatsFuns.jl/issues/115
-@testset "support of binomial distribution" begin
-    @test iszero(binompdf(1, 0.5, prevfloat(1.0)))
-    @test iszero(binompdf(1, 0.5, nextfloat(1.0)))
-    @test binomlogpdf(1, 0.5, prevfloat(1.0)) == -Inf
-    @test binomlogpdf(1, 0.5, nextfloat(1.0)) == -Inf
-end
+    # https://github.com/JuliaStats/StatsFuns.jl/issues/115
+    @testset "support of binomial distribution" begin
+        @test iszero(binompdf(1, 0.5, prevfloat(1.0)))
+        @test iszero(binompdf(1, 0.5, nextfloat(1.0)))
+        @test binomlogpdf(1, 0.5, prevfloat(1.0)) == -Inf
+        @test binomlogpdf(1, 0.5, nextfloat(1.0)) == -Inf
+    end
 
-@testset "binom special cases" begin
-    for (n, p, k) in ((5, 0.0, 0), (5, 1.0, 5))
-        @test iszero(binomlogpdf(n, p, k))
-        @test isone(binompdf(n, p, k))
+    @testset "binom special cases" begin
+        for (n, p, k) in ((5, 0.0, 0), (5, 1.0, 5))
+            @test iszero(binomlogpdf(n, p, k))
+            @test isone(binompdf(n, p, k))
+        end
     end
-end
 
-@testset "lstirling_asym" begin
-    # can test for equality here because the lhs is the way the value is created
-    @test Float32(lstirling_asym(1.0)) == @inferred lstirling_asym(1.0f0)
-    # for 64.0f0 the expansion is used but for 64.0 the BigFloat value is rounded
-    @test Float32(lstirling_asym(64.0)) ≈ @inferred lstirling_asym(64.0f0)
-end
+    @testset "lstirling_asym" begin
+        # can test for equality here because the lhs is the way the value is created
+        @test Float32(lstirling_asym(1.0)) == @inferred lstirling_asym(1.0f0)
+        # for 64.0f0 the expansion is used but for 64.0 the BigFloat value is rounded
+        @test Float32(lstirling_asym(64.0)) ≈ @inferred lstirling_asym(64.0f0)
+    end
 
-# https://github.com/JuliaStats/StatsFuns.jl/issues/143
-# https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/47
-@testset "betalogcdf: numerical issue" begin
-    # Mathematica: N[Log[CDF[BetaDistribution[6041, 2496], 1/10]], 10]
-    @test betalogcdf(6041, 2496, 0.1) ≈ -9020.029401
-    @test betainvlogcdf(6041, 2496, betalogcdf(6041, 2496, 0.1)) ≈ 0.1
-end
+    # https://github.com/JuliaStats/StatsFuns.jl/issues/143
+    # https://github.com/JuliaMath/HypergeometricFunctions.jl/issues/47
+    @testset "betalogcdf: numerical issue" begin
+        # Mathematica: N[Log[CDF[BetaDistribution[6041, 2496], 1/10]], 10]
+        @test betalogcdf(6041, 2496, 0.1) ≈ -9020.029401
+        @test betainvlogcdf(6041, 2496, betalogcdf(6041, 2496, 0.1)) ≈ 0.1
+    end
 
-# https://github.com/JuliaStats/StatsFuns.jl/issues/150
-@testset "gammalogcdf: numerical issue" begin
-    @test gammalogcdf(42648.50647826826, 2.2498007956420723e-5, 0.6991377135675367) ≈ -1933.269895904061741
-end
+    # https://github.com/JuliaStats/StatsFuns.jl/issues/150
+    @testset "gammalogcdf: numerical issue" begin
+        @test gammalogcdf(42648.50647826826, 2.2498007956420723e-5, 0.6991377135675367) ≈ -1933.269895904061741
+    end
 
-# https://github.com/JuliaStats/StatsFuns.jl/issues/154
-@testset "tvdistinvcdf: numerical issue" begin
-    @test isnan(@inferred(tdistinvcdf(0, 0.975)))
+    # https://github.com/JuliaStats/StatsFuns.jl/issues/154
+    @testset "tvdistinvcdf: numerical issue" begin
+        @test isnan(@inferred(tdistinvcdf(0, 0.975)))
+    end
 end
diff --git a/test/qa.jl b/test/qa.jl
index 85d2bcb..c9e7d46 100644
--- a/test/qa.jl
+++ b/test/qa.jl
@@ -1,10 +1,10 @@
-using StatsFuns
-using Aqua: Aqua
-using ExplicitImports: ExplicitImports
-using JET: JET
-using Test
+@testitem "QA" begin
+    using StatsFuns
+    using Aqua: Aqua
+    using ExplicitImports: ExplicitImports
+    using JET: JET
+    using Test
 
-@testset "QA" begin
     @testset "Aqua" begin
         Aqua.test_all(StatsFuns)
     end
diff --git a/test/rmath.jl b/test/rmath.jl
index 8d5bd60..9cea300 100644
--- a/test/rmath.jl
+++ b/test/rmath.jl
@@ -1,153 +1,160 @@
-using StatsFuns
-using Rmath: Rmath
-using Test
-
-include("utils.jl")
-
-function check_rmath(statsfun, rmathfun, params, a, isprob, rtol)
-    v = @inferred(statsfun(params..., a))
-    rv = @inferred(rmathfun(a, params...))
-    @test v isa float(Base.promote_typeof(params..., a))
-    return if isprob
-        @test v ≈ oftype(v, rv) rtol = rtol nans = true
-    else
-        @test v ≈ oftype(v, rv) atol = rtol rtol = rtol nans = true
-    end
-end
+@testsnippet RMathHelpers begin
+    using Rmath: Rmath
 
-function rmathcomp(basename::String, params, X::AbstractArray, rtol = _default_rtol(params, X))
-    rbasename = if basename == "nfdist"
-        "nf"
-    elseif basename == "ntdist"
-        "nt"
-    elseif basename == "srdist"
-        "tukey"
-    else
-        basename
+    # default relative tolerance for comparisons with Rmath
+    function _default_rtol(params, X::AbstractArray)
+        return _default_rtol(float(promote_type(Base.promote_typeof(params...), eltype(X))))
     end
-
-    if isdefined(Rmath, Symbol(:d, rbasename))
-        stats_pdf = getproperty(@__MODULE__, Symbol(basename, :pdf))
-        rmath_pdf = let f = getproperty(Rmath, Symbol(:d, rbasename))
-            (a, params...) -> f(a, params..., false)
-        end
-        @testset "pdf with x=$x" for x in X
-            check_rmath(
-                stats_pdf, rmath_pdf,
-                params, x, true, rtol
-            )
-        end
-
-        stats_logpdf = getproperty(@__MODULE__, Symbol(basename, :logpdf))
-        rmath_logpdf = let f = getproperty(Rmath, Symbol(:d, rbasename))
-            (a, params...) -> f(a, params..., true)
-        end
-        @testset "logpdf with x=$x" for x in X
-            check_rmath(
-                stats_logpdf, rmath_logpdf,
-                params, x, false, rtol
-            )
+    _default_rtol(::Type{Float64}) = eps(Float64)^(7 // 8)
+    _default_rtol(::Type{Float32}) = eps(Float32)^(3 // 4)
+    _default_rtol(::Type{Float16}) = eps(Float16)^(2 // 3)
+
+    function check_rmath(statsfun, rmathfun, params, a, isprob, rtol)
+        v = @inferred(statsfun(params..., a))
+        rv = @inferred(rmathfun(a, params...))
+        @test v isa float(Base.promote_typeof(params..., a))
+        return if isprob
+            @test v ≈ oftype(v, rv) rtol = rtol nans = true
+        else
+            @test v ≈ oftype(v, rv) atol = rtol rtol = rtol nans = true
         end
     end
 
-    if isdefined(Rmath, Symbol(:p, rbasename))
-        stats_cdf = getproperty(@__MODULE__, Symbol(basename, :cdf))
-        rmath_cdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
-            (a, params...) -> f(a, params..., true, false)
-        end
-        @testset "cdf with x=$x" for x in X
-            check_rmath(stats_cdf, rmath_cdf, params, x, true, rtol)
-        end
-
-        stats_ccdf = getproperty(@__MODULE__, Symbol(basename, :ccdf))
-        rmath_ccdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
-            (a, params...) -> f(a, params..., false, false)
-        end
-        @testset "ccdf with x=$x" for x in X
-            check_rmath(stats_ccdf, rmath_ccdf, params, x, true, rtol)
+    function rmathcomp(basename::String, params, X::AbstractArray, rtol = _default_rtol(params, X))
+        rbasename = if basename == "nfdist"
+            "nf"
+        elseif basename == "ntdist"
+            "nt"
+        elseif basename == "srdist"
+            "tukey"
+        else
+            basename
         end
 
-        stats_logcdf = getproperty(@__MODULE__, Symbol(basename, :logcdf))
-        rmath_logcdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
-            (a, params...) -> f(a, params..., true, true)
-        end
-        @testset "logcdf with x=$x" for x in X
-            check_rmath(stats_logcdf, rmath_logcdf, params, x, false, rtol)
-        end
+        if isdefined(Rmath, Symbol(:d, rbasename))
+            stats_pdf = getproperty(@__MODULE__, Symbol(basename, :pdf))
+            rmath_pdf = let f = getproperty(Rmath, Symbol(:d, rbasename))
+                (a, params...) -> f(a, params..., false)
+            end
+            @testset "pdf with x=$x" for x in X
+                check_rmath(
+                    stats_pdf, rmath_pdf,
+                    params, x, true, rtol
+                )
+            end
 
-        stats_logccdf = getproperty(@__MODULE__, Symbol(basename, :logccdf))
-        rmath_logccdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
-            (a, params...) -> f(a, params..., false, true)
-        end
-        @testset "logccdf with x=$x" for x in X
-            check_rmath(stats_logccdf, rmath_logccdf, params, x, false, rtol)
+            stats_logpdf = getproperty(@__MODULE__, Symbol(basename, :logpdf))
+            rmath_logpdf = let f = getproperty(Rmath, Symbol(:d, rbasename))
+                (a, params...) -> f(a, params..., true)
+            end
+            @testset "logpdf with x=$x" for x in X
+                check_rmath(
+                    stats_logpdf, rmath_logpdf,
+                    params, x, false, rtol
+                )
+            end
         end
 
-        #=
-        R version of signrank has system variation
-        julia> psignrank(18,10,false,true) # windows
-        -0.2076393647782445
-        julia> psignrank(18,10,false,true) # linux
-        -0.20763936477824452
-        This slight difference causes test failures for the inverse functions,
-        due to a slight shift in the location of the discontinuity.
-    
-        This also holds true for wilcox.
-        =#
-        test_inv = (basename != "signrank" && basename != "wilcox") || !Sys.islinux()
-        if isdefined(Rmath, Symbol(:q, rbasename)) && test_inv
-            stats_invcdf = getproperty(@__MODULE__, Symbol(basename, :invcdf))
-            rmath_invcdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+        if isdefined(Rmath, Symbol(:p, rbasename))
+            stats_cdf = getproperty(@__MODULE__, Symbol(basename, :cdf))
+            rmath_cdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
                 (a, params...) -> f(a, params..., true, false)
             end
-            p = rmath_cdf.(X, params...)
-            @testset "invcdf with q=$_p" for _p in p
-                check_rmath(stats_invcdf, rmath_invcdf, params, _p, false, rtol)
+            @testset "cdf with x=$x" for x in X
+                check_rmath(stats_cdf, rmath_cdf, params, x, true, rtol)
             end
 
-            stats_invccdf = getproperty(@__MODULE__, Symbol(basename, :invccdf))
-            rmath_invccdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+            stats_ccdf = getproperty(@__MODULE__, Symbol(basename, :ccdf))
+            rmath_ccdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
                 (a, params...) -> f(a, params..., false, false)
             end
-            cp = rmath_ccdf.(X, params...)
-            @testset "invccdf with q=$_p" for _p in cp
-                check_rmath(stats_invccdf, rmath_invccdf, params, _p, false, rtol)
+            @testset "ccdf with x=$x" for x in X
+                check_rmath(stats_ccdf, rmath_ccdf, params, x, true, rtol)
             end
 
-            stats_invlogcdf = getproperty(@__MODULE__, Symbol(basename, :invlogcdf))
-            rmath_invlogcdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+            stats_logcdf = getproperty(@__MODULE__, Symbol(basename, :logcdf))
+            rmath_logcdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
                 (a, params...) -> f(a, params..., true, true)
             end
-            lp = rmath_logcdf.(X, params...)
-            @testset "invlogcdf with log(q)=$_p" for _p in lp
-                check_rmath(stats_invlogcdf, rmath_invlogcdf, params, _p, false, rtol)
+            @testset "logcdf with x=$x" for x in X
+                check_rmath(stats_logcdf, rmath_logcdf, params, x, false, rtol)
             end
 
-            stats_invlogccdf = getproperty(@__MODULE__, Symbol(basename, :invlogccdf))
-            rmath_invlogccdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+            stats_logccdf = getproperty(@__MODULE__, Symbol(basename, :logccdf))
+            rmath_logccdf = let f = getproperty(Rmath, Symbol(:p, rbasename))
                 (a, params...) -> f(a, params..., false, true)
             end
-            lcp = rmath_logccdf.(X, params...)
-            @testset "invlogccdf with log(q)=$_p" for _p in lcp
-                check_rmath(stats_invlogccdf, rmath_invlogccdf, params, _p, false, rtol)
+            @testset "logccdf with x=$x" for x in X
+                check_rmath(stats_logccdf, rmath_logccdf, params, x, false, rtol)
+            end
+
+            #=
+            R version of signrank has system variation
+            julia> psignrank(18,10,false,true) # windows
+            -0.2076393647782445
+            julia> psignrank(18,10,false,true) # linux
+            -0.20763936477824452
+            This slight difference causes test failures for the inverse functions,
+            due to a slight shift in the location of the discontinuity.
+
+            This also holds true for wilcox.
+            =#
+            test_inv = (basename != "signrank" && basename != "wilcox") || !Sys.islinux()
+            if isdefined(Rmath, Symbol(:q, rbasename)) && test_inv
+                stats_invcdf = getproperty(@__MODULE__, Symbol(basename, :invcdf))
+                rmath_invcdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+                    (a, params...) -> f(a, params..., true, false)
+                end
+                p = rmath_cdf.(X, params...)
+                @testset "invcdf with q=$_p" for _p in p
+                    check_rmath(stats_invcdf, rmath_invcdf, params, _p, false, rtol)
+                end
+
+                stats_invccdf = getproperty(@__MODULE__, Symbol(basename, :invccdf))
+                rmath_invccdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+                    (a, params...) -> f(a, params..., false, false)
+                end
+                cp = rmath_ccdf.(X, params...)
+                @testset "invccdf with q=$_p" for _p in cp
+                    check_rmath(stats_invccdf, rmath_invccdf, params, _p, false, rtol)
+                end
+
+                stats_invlogcdf = getproperty(@__MODULE__, Symbol(basename, :invlogcdf))
+                rmath_invlogcdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+                    (a, params...) -> f(a, params..., true, true)
+                end
+                lp = rmath_logcdf.(X, params...)
+                @testset "invlogcdf with log(q)=$_p" for _p in lp
+                    check_rmath(stats_invlogcdf, rmath_invlogcdf, params, _p, false, rtol)
+                end
+
+                stats_invlogccdf = getproperty(@__MODULE__, Symbol(basename, :invlogccdf))
+                rmath_invlogccdf = let f = getproperty(Rmath, Symbol(:q, rbasename))
+                    (a, params...) -> f(a, params..., false, true)
+                end
+                lcp = rmath_logccdf.(X, params...)
+                @testset "invlogccdf with log(q)=$_p" for _p in lcp
+                    check_rmath(stats_invlogccdf, rmath_invlogccdf, params, _p, false, rtol)
+                end
             end
         end
-    end
 
-    return nothing
-end
+        return nothing
+    end
 
-function rmathcomp_tests(basename::String, configs)
-    return @testset "$basename" begin
-        @testset "params: $params" for (params, data) in configs
-            rmathcomp(basename, params, data)
+    function rmathcomp_tests(basename::String, configs)
+        return @testset "$basename" begin
+            @testset "params: $params" for (params, data) in configs
+                rmathcomp(basename, params, data)
+            end
         end
     end
 end
 
-### Test cases
+@testitem "RMath beta" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
-@testset "RMath" begin
     rmathcomp_tests(
         "beta", [
             ((0.1, 1.0), 0.0:0.01:1.0),
@@ -219,6 +226,20 @@ end
         end
     end
 
+    # Test values outside of the support
+    rmathcomp_tests(
+        "beta", [
+            ((1.0, 1.0), [-10.0, -6.3, 2.1, 23.5]),
+            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10, 21 // 10, 47 // 2]),
+            ((1, 1), [-10, -6, 2, 24]),
+        ]
+    )
+end
+
+@testitem "RMath binom" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "binom", [
             ((1, 0.5), 0.0:1.0),
@@ -236,6 +257,20 @@ end
         ]
     )
 
+    # Test values outside of the support
+    rmathcomp_tests(
+        "binom", [
+            ((5, 0.5), [-8, -2.3, 1.2, 5.4, 11.9]),
+            ((5, 1 // 2), [-8, -23 // 10, 6 // 5, 27 // 5, 119 // 10]),
+            ((5, 1 // 2), [-8, -2, 6, 12]),
+        ]
+    )
+end
+
+@testitem "RMath chisq" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "chisq", [
             ((1,), 0.0:0.1:8.0),
@@ -248,6 +283,11 @@ end
             ((1), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath fdist" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "fdist", [
@@ -263,6 +303,20 @@ end
         ]
     )
 
+    # Test values outside of the support
+    rmathcomp_tests(
+        "fdist", [
+            ((1.0, 1.0), [-10.0, -6.3]),
+            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10]),
+            ((1, 1), [-10, -6]),
+        ]
+    )
+end
+
+@testitem "RMath gamma" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "gamma", [
             ((1.0, 1.0), (0.0:0.05:12.0)),
@@ -281,6 +335,20 @@ end
         ]
     )
 
+    # Test values outside of the support
+    rmathcomp_tests(
+        "gamma", [
+            ((1.0, 1.0), [-10.0, -6.3]),
+            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10]),
+            ((1, 1), [-10, -6]),
+        ]
+    )
+end
+
+@testitem "RMath hyper" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "hyper", [
             ((2, 3, 4), 0.0:4.0),
@@ -290,6 +358,11 @@ end
             ((2, 3, 4), (0 // 1):(4 // 1)),
         ]
     )
+end
+
+@testitem "RMath nbeta" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "nbeta", [
@@ -306,6 +379,11 @@ end
             ((1.0, 1.0, 0.0), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath nbinom" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "nbinom", [
@@ -322,6 +400,11 @@ end
             ((1, 0.5f0), [-Inf32, Inf32, NaN32, -2, 1.1]),
         ]
     )
+end
+
+@testitem "RMath nchisq" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "nchisq", [
@@ -336,6 +419,11 @@ end
             ((2, 1), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath nfdist" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "nfdist", [
@@ -351,6 +439,11 @@ end
             ((1.0, 1.0, 0.0), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath norm" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "norm", [
@@ -371,6 +464,11 @@ end
             #((0f0, 1f0), -Float16(6):Float16(0.01):Float16(6)),
         ]
     )
+end
+
+@testitem "RMath ntdist" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "ntdist", [
@@ -385,6 +483,11 @@ end
             ((10, 1), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath pois" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "pois", [
@@ -401,6 +504,20 @@ end
         ]
     )
 
+    # Test values outside of the support
+    rmathcomp_tests(
+        "pois", [
+            ((0.5,), [-10, -2.5, 1.3, 8.7]),
+            ((1 // 2,), [-10, -5 // 2, 13 // 10, 87 // 10]),
+            ((1,), [-10, -3]),
+        ]
+    )
+end
+
+@testitem "RMath tdist" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "tdist", [
             ((1,), -5.0:0.1:5.0),
@@ -415,6 +532,11 @@ end
             ((1,), [-Inf, Inf]),
         ]
     )
+end
+
+@testitem "RMath signrank" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "signrank", [
@@ -433,6 +555,11 @@ end
     @test signrankinvlogccdf.(50, signranklogccdf.(50, -1:1274)) == [0; 0:1274]
     @test isnan(signrankinvlogccdf.(50, signranklogccdf.(50, 1275)))
     @test isnan(signrankinvlogccdf.(50, signranklogccdf.(50, 1276)))
+end
+
+@testitem "RMath srdist" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
 
     rmathcomp_tests(
         "srdist", [
@@ -447,6 +574,19 @@ end
         ]
     )
 
+    # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5.
+    # Reduced df or k allows convergence. This test documents this behavior.
+    x = 0.15
+    q = srdistcdf(10, 5, 0.15)
+    rx = srdistinvcdf(10, 5, q)
+    rtol = 100eps(1.0)
+    @test_broken x ≈ rx atol = rtol rtol = rtol nans = true
+end
+
+@testitem "RMath wilcox" setup=[RMathHelpers] begin
+    using StatsFuns
+    using Test
+
     rmathcomp_tests(
         "wilcox", [
             ((3, 4), -2:13),
@@ -458,7 +598,6 @@ end
         ]
     )
 
-
     @test wilcoxinvcdf.(10, 10, wilcoxcdf.(10, 10, -1:101)) == [0; 0:100; 100]
     @test wilcoxinvccdf.(10, 10, wilcoxccdf.(10, 10, -1:101)) == [0; 0:100; 100]
     @test wilcoxinvlogcdf.(10, 10, wilcoxlogcdf.(10, 10, 0:101)) == [0:100; 100]
@@ -466,49 +605,4 @@ end
     @test wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, -1:99)) == [0; 0:99]
     @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 100)))
     @test isnan(wilcoxinvlogccdf.(10, 10, wilcoxlogccdf.(10, 10, 101)))
-
-    # Note: Convergence fails in srdist with cdf values below 0.16 with df = 10, k = 5.
-    # Reduced df or k allows convergence. This test documents this behavior.
-    x = 0.15
-    q = srdistcdf(10, 5, 0.15)
-    rx = srdistinvcdf(10, 5, q)
-    rtol = 100eps(1.0)
-    @test_broken x ≈ rx atol = rtol rtol = rtol nans = true
-
-    # Test values outside of the support
-    rmathcomp_tests(
-        "beta", [
-            ((1.0, 1.0), [-10.0, -6.3, 2.1, 23.5]),
-            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10, 21 // 10, 47 // 2]),
-            ((1, 1), [-10, -6, 2, 24]),
-        ]
-    )
-    rmathcomp_tests(
-        "binom", [
-            ((5, 0.5), [-8, -2.3, 1.2, 5.4, 11.9]),
-            ((5, 1 // 2), [-8, -23 // 10, 6 // 5, 27 // 5, 119 // 10]),
-            ((5, 1 // 2), [-8, -2, 6, 12]),
-        ]
-    )
-    rmathcomp_tests(
-        "fdist", [
-            ((1.0, 1.0), [-10.0, -6.3]),
-            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10]),
-            ((1, 1), [-10, -6]),
-        ]
-    )
-    rmathcomp_tests(
-        "gamma", [
-            ((1.0, 1.0), [-10.0, -6.3]),
-            ((1 // 1, 1 // 1), [-10 // 1, -63 // 10]),
-            ((1, 1), [-10, -6]),
-        ]
-    )
-    rmathcomp_tests(
-        "pois", [
-            ((0.5,), [-10, -2.5, 1.3, 8.7]),
-            ((1 // 2,), [-10, -5 // 2, 13 // 10, 87 // 10]),
-            ((1,), [-10, -3]),
-        ]
-    )
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0cd38a1..2380de7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,2 @@
-tests = ["rmath", "generic", "misc", "chainrules", "inverse", "tvpack", "qa"]
-
-for t in tests
-    fp = "$t.jl"
-    println("* running $fp")
-    include(fp)
-end
+using TestItemRunner
+@run_package_tests
diff --git a/test/tvpack.jl b/test/tvpack.jl
index a25ae42..6dda83e 100644
--- a/test/tvpack.jl
+++ b/test/tvpack.jl
@@ -1,30 +1,32 @@
-using StatsFuns, Test
+@testitem "TVPack" begin
+    using StatsFuns, Test
 
-@testset "bvncdf: bivariate normal cdf" begin
+    @testset "bvncdf: bivariate normal cdf" begin
 
-    @test bvncdf(0.0, 0.0, 0.0) == 0.25
-    @test bvncdf(0.0, Inf, 0.0) == 0.5
-    @test bvncdf(Inf, 0.0, 0.0) == 0.5
-    @test bvncdf(Inf, Inf, 0.0) == 1.0
+        @test bvncdf(0.0, 0.0, 0.0) == 0.25
+        @test bvncdf(0.0, Inf, 0.0) == 0.5
+        @test bvncdf(Inf, 0.0, 0.0) == 0.5
+        @test bvncdf(Inf, Inf, 0.0) == 1.0
 
-    for r in -1.0:0.25:1.0
-        for x in -10.0:0.5:10.0
-            @test bvncdf(x, Inf, r) == normcdf(x)
-            @test bvncdf(Inf, x, r) == normcdf(x)
-            @test bvncdf(x, -Inf, r) == 0.0
-            @test bvncdf(-Inf, x, r) == 0.0
+        for r in -1.0:0.25:1.0
+            for x in -10.0:0.5:10.0
+                @test bvncdf(x, Inf, r) == normcdf(x)
+                @test bvncdf(Inf, x, r) == normcdf(x)
+                @test bvncdf(x, -Inf, r) == 0.0
+                @test bvncdf(-Inf, x, r) == 0.0
+            end
         end
-    end
 
-    @test_throws DomainError bvncdf(0.0, 0.0, -2.0)
-    @test_throws DomainError bvncdf(0.0, 0.0, 2.0)
+        @test_throws DomainError bvncdf(0.0, 0.0, -2.0)
+        @test_throws DomainError bvncdf(0.0, 0.0, 2.0)
 
-    @test bvncdf(0.0, -100000.0, 0.0) ≈ 0.0
-    @test bvncdf(0.0, 100000.0, 0.0) ≈ 0.5
-    @test bvncdf(100000.0, 100000.0, 0.0) ≈ 1.0
-    @test bvncdf(-100000.0, 100000.0, 0.0) ≈ 0.0
+        @test bvncdf(0.0, -100000.0, 0.0) ≈ 0.0
+        @test bvncdf(0.0, 100000.0, 0.0) ≈ 0.5
+        @test bvncdf(100000.0, 100000.0, 0.0) ≈ 1.0
+        @test bvncdf(-100000.0, 100000.0, 0.0) ≈ 0.0
 
-    for x in -100.0:100.0, y in -100.0:100.0, r in -1.0:0.1:1.0
-        @test !isnan(bvncdf(x, y, r))
+        for x in -100.0:100.0, y in -100.0:100.0, r in -1.0:0.1:1.0
+            @test !isnan(bvncdf(x, y, r))
+        end
     end
 end
diff --git a/test/utils.jl b/test/utils.jl
index 19ba4d9..2cdb891 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1,3 +1,5 @@
+@testmodule Utils begin
+
 # default relative tolerance for comparisons with Rmath
 function _default_rtol(params, X::AbstractArray)
     # has to take into account `params` as well since otherwise e.g. `X::Array{<:Rational}`
@@ -12,3 +14,5 @@ end
 _default_rtol(::Type{Float64}) = eps(Float64)^(7 // 8) # ~2.0e-14
 _default_rtol(::Type{Float32}) = eps(Float32)^(3 // 4) # ~6.4e-6
 _default_rtol(::Type{Float16}) = eps(Float16)^(2 // 3) # ~9.9e-3
+
+end