diff --git a/Project.toml b/Project.toml
index e707ba77..fcdef390 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,22 +1,25 @@
 name = "HypothesisTests"
 uuid = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
-version = "0.10.13"
+version = "0.11.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Combinatorics = "1"
-Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
+Distributions = "0.25.90"
 Rmath = "0.5, 0.6, 0.7"
 Roots = "0.7, 0.8, 1.0, 2"
+StatsAPI = "1.6"
 StatsBase = "0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34"
 julia = "1.3"
 
diff --git a/src/HypothesisTests.jl b/src/HypothesisTests.jl
index 0356237a..e3d12e2b 100644
--- a/src/HypothesisTests.jl
+++ b/src/HypothesisTests.jl
@@ -28,73 +28,17 @@ using Statistics, Random, LinearAlgebra
 using Distributions, Roots, StatsBase
 using Combinatorics: combinations, permutations
 using Rmath: pwilcox, psignrank
+using Printf: @printf
 
-import StatsBase.confint
+import StatsAPI
+using StatsAPI: HypothesisTest, confint, pvalue
 
-export testname, pvalue, confint
-abstract type HypothesisTest end
+export testname, pvalue, confint, dof, nobs
 
 check_same_length(x::AbstractVector, y::AbstractVector) = if length(x) != length(y)
     throw(DimensionMismatch("Vectors must be the same length"))
 end
 
-"""
-    confint(test::HypothesisTest; level = 0.95, tail = :both)
-
-Compute a confidence interval C with coverage `level`.
-
-If `tail` is `:both` (default), then a two-sided confidence interval is returned. If `tail`
-is `:left` or `:right`, then a one-sided confidence interval is returned.
-
-!!! note
-    Most of the implemented confidence intervals are *strongly consistent*, that is, the
-    confidence interval with coverage `level` does not contain the test statistic under
-    ``h_0`` if and only if the corresponding test rejects the null hypothesis
-    ``h_0: θ = θ_0``:
-    ```math
-        C (x, level) = \\{θ : p_θ (x) > 1 - level\\},
-    ```
-    where ``p_θ`` is the [`pvalue`](@ref) of the corresponding test.
-"""
-function confint end
-
-"""
-    pvalue(test::HypothesisTest; tail = :both)
-
-Compute the p-value for a given significance test.
-
-If `tail` is `:both` (default), then the p-value for the two-sided test is returned. If
-`tail` is `:left` or `:right`, then a one-sided test is performed.
-"""
-function pvalue end
-
-# Basic function for finding a p-value given a distribution and tail
-function pvalue(dist::ContinuousUnivariateDistribution, x::Number; tail=:both)
-    check_tail(tail)
-
-    if tail == :both
-        p = 2 * min(cdf(dist, x), ccdf(dist, x))
-        min(p, oneunit(p)) # if P(X = x) > 0, then possibly p > 1
-    elseif tail == :left
-        cdf(dist, x)
-    else # tail == :right
-        ccdf(dist, x)
-    end
-end
-
-function pvalue(dist::DiscreteUnivariateDistribution, x::Number; tail=:both)
-    check_tail(tail)
-
-    if tail == :both
-        p = 2 * min(ccdf(dist, x-1), cdf(dist, x))
-        min(p, oneunit(p)) # if P(X = x) > 0, then possibly p > 1
-    elseif tail == :left
-        cdf(dist, x)
-    else # tail == :right
-        ccdf(dist, x-1)
-    end
-end
-
 function check_level(level::Float64)
     if level >= 1 || level <= 0.5
         throw(ArgumentError("coverage level $level not in range (0.5, 1)"))
@@ -114,7 +58,7 @@ function Base.show(_io::IO, test::T) where T<:HypothesisTest
     println(io, repeat("-", length(testname(test))))
 
     # population details
-    has_ci = applicable(StatsBase.confint, test)
+    has_ci = applicable(confint, test)
     (param_name, param_under_h0, param_estimate) = population_param_of_interest(test)
     println(io, "Population details:")
     println(io, "    parameter of interest:   $param_name")
@@ -126,7 +70,7 @@ function Base.show(_io::IO, test::T) where T<:HypothesisTest
     println(io)
 
     if has_ci
-        ci = map(x -> round.(x; sigdigits=4, base=10), StatsBase.confint(test))
+        ci = map(x -> round.(x; sigdigits=4, base=10), confint(test))
         print(io, "    95% confidence interval: ")
         show(io, ci)
         println(io)
@@ -176,7 +120,6 @@ function show_params(io::IO, test::T, ident="") where T<:HypothesisTest
     end
 end
 
-include("deprecated.jl")
 include("common.jl")
 
 include("binomial.jl")
diff --git a/src/anderson_darling.jl b/src/anderson_darling.jl
index d2e26688..807f396d 100644
--- a/src/anderson_darling.jl
+++ b/src/anderson_darling.jl
@@ -52,7 +52,7 @@ function show_params(io::IO, x::OneSampleADTest, ident = "")
 end
 
 ### G. and J. Marsaglia, "Evaluating the Anderson-Darling Distribution", Journal of Statistical Software, 2004
-function pvalue(t::OneSampleADTest)
+function StatsAPI.pvalue(t::OneSampleADTest)
     g1(x) = sqrt(x)*(1.0-x)*(49.0x-102.0)
     g2(x) = -0.00022633 + (6.54034 - (14.6538 - (14.458 - (8.259 - 1.91864x)x)x)x)x
     g3(x) = -130.2137 + (745.2337 - (1705.091 - (1950.646 - (1116.360 - 255.7844x)x)x)x)x
@@ -238,7 +238,7 @@ function pvalueasym(x::KSampleADTest)
     return exp(lp0)/(1 + exp(lp0))
 end
 
-pvalue(x::KSampleADTest) = x.nsim == 0 ? pvalueasym(x) : pvaluesim(x)
+StatsAPI.pvalue(x::KSampleADTest) = x.nsim == 0 ? pvalueasym(x) : pvaluesim(x)
 
 function adkvals(Z⁺, N, samples)
     k = length(samples)
diff --git a/src/augmented_dickey_fuller.jl b/src/augmented_dickey_fuller.jl
index 67e45e30..27f950da 100644
--- a/src/augmented_dickey_fuller.jl
+++ b/src/augmented_dickey_fuller.jl
@@ -217,4 +217,4 @@ function show_params(io::IO, x::ADFTest, ident)
     println(io)
 end
 
-pvalue(x::ADFTest) = HypothesisTests.pvalue(Normal(0, 1), adf_pv_aux(x.stat, x.deterministic); tail=:left)
+StatsAPI.pvalue(x::ADFTest) = pvalue(Normal(0, 1), adf_pv_aux(x.stat, x.deterministic); tail=:left)
diff --git a/src/bartlett.jl b/src/bartlett.jl
index 944f70d8..263ea0fd 100644
--- a/src/bartlett.jl
+++ b/src/bartlett.jl
@@ -45,12 +45,12 @@ function BartlettTest(X::AbstractMatrix, Y::AbstractMatrix)
     return BartlettTest(L′, p, nx, ny)
 end
 
-StatsBase.nobs(B::BartlettTest) = (B.nx, B.ny)
-StatsBase.dof(B::BartlettTest) = div(B.p * (B.p + 1), 2)
+StatsAPI.nobs(B::BartlettTest) = (B.nx, B.ny)
+StatsAPI.dof(B::BartlettTest) = div(B.p * (B.p + 1), 2)
 
 testname(::BartlettTest) = "Bartlett's Test for Equality of Covariance Matrices"
 default_tail(::BartlettTest) = :right
-pvalue(B::BartlettTest; tail=:right) = pvalue(Chisq(dof(B)), B.L′, tail=tail)
+StatsAPI.pvalue(B::BartlettTest; tail=:right) = pvalue(Chisq(dof(B)), B.L′, tail=tail)
 
 function show_params(io::IO, B::BartlettTest, indent="")
     println(io, indent, "number of observations: ", nobs(B))
diff --git a/src/binomial.jl b/src/binomial.jl
index d2aa74f4..2e3be219 100644
--- a/src/binomial.jl
+++ b/src/binomial.jl
@@ -66,7 +66,7 @@ function show_params(io::IO, x::BinomialTest, ident="")
     println(io, ident, "number of successes:    $(x.x)")
 end
 
-pvalue(x::BinomialTest; tail=:both) = pvalue(Binomial(x.n, x.p), x.x; tail=tail)
+StatsAPI.pvalue(x::BinomialTest; tail=:both) = pvalue(Binomial(x.n, x.p), x.x; tail=tail)
 
 # Confidence interval
 """
@@ -102,13 +102,13 @@ of the following methods. Possible values for `method` are:
   * [Binomial confidence interval on Wikipedia](https://en.wikipedia.org/wiki/
     Binomial_proportion_confidence_interval)
 """
-function StatsBase.confint(x::BinomialTest; level::Float64=0.95, tail=:both, method=:clopper_pearson)
+function StatsAPI.confint(x::BinomialTest; level::Float64=0.95, tail=:both, method=:clopper_pearson)
     check_level(level)
 
     if tail == :left
-        (0.0, StatsBase.confint(x, level=1-(1-level)*2, method=method)[2])
+        (0.0, confint(x, level=1-(1-level)*2, method=method)[2])
     elseif tail == :right
-        (StatsBase.confint(x, level=1-(1-level)*2, method=method)[1], 1.0)
+        (confint(x, level=1-(1-level)*2, method=method)[1], 1.0)
     elseif tail == :both
         if method == :clopper_pearson
             ci_clopper_pearson(x, 1-level)
@@ -218,9 +218,9 @@ function show_params(io::IO, x::SignTest, ident="")
     println(io, ident, text2, x.x)
 end
 
-pvalue(x::SignTest; tail=:both) = pvalue(Binomial(x.n, 0.5), x.x; tail=tail)
+StatsAPI.pvalue(x::SignTest; tail=:both) = pvalue(Binomial(x.n, 0.5), x.x; tail=tail)
 
-function StatsBase.confint(x::SignTest; level::Float64=0.95, tail=:both)
+function StatsAPI.confint(x::SignTest; level::Float64=0.95, tail=:both)
     check_level(level)
 
     if tail == :left
diff --git a/src/box_test.jl b/src/box_test.jl
index 2d4a3e9a..6e06e933 100644
--- a/src/box_test.jl
+++ b/src/box_test.jl
@@ -71,7 +71,7 @@ function show_params(io::IO, x::BoxPierceTest, ident)
     println(io, ident, "Q statistic:                    ", x.Q)
 end
 
-pvalue(x::BoxPierceTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
+StatsAPI.pvalue(x::BoxPierceTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
 
 #Ljung-Box test
 
@@ -120,4 +120,4 @@ function show_params(io::IO, x::LjungBoxTest, ident)
     println(io, ident, "Q statistic:                    ", x.Q)
 end
 
-pvalue(x::LjungBoxTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
+StatsAPI.pvalue(x::LjungBoxTest) = pvalue(Chisq(x.lag-x.dof), x.Q; tail=:right)
diff --git a/src/breusch_godfrey.jl b/src/breusch_godfrey.jl
index 21be0d03..d5b5b767 100644
--- a/src/breusch_godfrey.jl
+++ b/src/breusch_godfrey.jl
@@ -76,4 +76,4 @@ function show_params(io::IO, x::BreuschGodfreyTest, ident)
     println(io, ident, "T*R^2 statistic:                ", x.BG)
 end
 
-pvalue(x::BreuschGodfreyTest) = pvalue(Chisq(x.lag), x.BG; tail=:right)
+StatsAPI.pvalue(x::BreuschGodfreyTest) = pvalue(Chisq(x.lag), x.BG; tail=:right)
diff --git a/src/circular.jl b/src/circular.jl
index 248faa5f..73d509ad 100644
--- a/src/circular.jl
+++ b/src/circular.jl
@@ -60,7 +60,7 @@ function show_params(io::IO, x::RayleighTest, ident="")
     println(io, ident, "test statistic:         $(x.Rbar^2 * x.n)")
 end
 
-function pvalue(x::RayleighTest)
+function StatsAPI.pvalue(x::RayleighTest)
     Z = x.Rbar^2 * x.n
     x.n > 1e6 ? exp(-Z) :
         exp(-Z)*(1+(2*Z-Z^2)/(4*x.n)-(24*Z - 132*Z^2 + 76*Z^3 - 9*Z^4)/(288*x.n^2))
@@ -135,7 +135,7 @@ function tlinear_Z(x::FisherTLinearAssociation)
 end
 
 # p-values
-function pvalue(x::FisherTLinearAssociation; tail=:both)
+function StatsAPI.pvalue(x::FisherTLinearAssociation; tail=:both)
     n = length(x.theta)
     if n == 0
         return NaN
@@ -218,7 +218,7 @@ function show_params(io::IO, x::JammalamadakaCircularCorrelation, ident="")
     println(io, ident, "test statistic: $(x.Z)")
 end
 
-pvalue(x::JammalamadakaCircularCorrelation; tail=:both) = pvalue(Normal(), x.Z; tail=tail)
+StatsAPI.pvalue(x::JammalamadakaCircularCorrelation; tail=:both) = pvalue(Normal(), x.Z; tail=tail)
 
 ## GENERAL
 
diff --git a/src/correlation.jl b/src/correlation.jl
index 4bbf3233..40db405a 100644
--- a/src/correlation.jl
+++ b/src/correlation.jl
@@ -67,10 +67,10 @@ function population_param_of_interest(p::CorrelationTest)
     (param, zero(p.r), p.r)
 end
 
-StatsBase.nobs(p::CorrelationTest) = p.n
-StatsBase.dof(p::CorrelationTest) = p.n - 2 - p.k
+StatsAPI.nobs(p::CorrelationTest) = p.n
+StatsAPI.dof(p::CorrelationTest) = p.n - 2 - p.k
 
-function StatsBase.confint(test::CorrelationTest{T}, level::Float64=0.95) where T
+function StatsAPI.confint(test::CorrelationTest{T}, level::Float64=0.95) where T
     dof(test) > 1 || return (-one(T), one(T))  # Otherwise we can get NaNs
     q = quantile(Normal(), 1 - (1-level) / 2)
     fisher = atanh(test.r)
@@ -81,7 +81,7 @@ function StatsBase.confint(test::CorrelationTest{T}, level::Float64=0.95) where
 end
 
 default_tail(::CorrelationTest) = :both
-pvalue(test::CorrelationTest; tail=:both) = pvalue(TDist(dof(test)), test.t, tail=tail)
+StatsAPI.pvalue(test::CorrelationTest; tail=:both) = pvalue(TDist(dof(test)), test.t, tail=tail)
 
 function show_params(io::IO, test::CorrelationTest, indent="")
     println(io, indent, "number of observations:          ", nobs(test))
diff --git a/src/deprecated.jl b/src/deprecated.jl
deleted file mode 100644
index b88a55b9..00000000
--- a/src/deprecated.jl
+++ /dev/null
@@ -1,12 +0,0 @@
-using Base: @deprecate
-
-@deprecate ci(args...) confint(args...)
-
-@deprecate MultinomialLRT MultinomialLRTest
-@deprecate OneSampleHotellingT2 OneSampleHotellingT2Test
-@deprecate EqualCovHotellingT2 EqualCovHotellingT2Test
-@deprecate UnequalCovHotellingT2 UnequalCovHotellingT2Test
-
-@deprecate BartlettsTest BartlettTest
-
-@deprecate confint(x::HypothesisTest, alpha::Real; kwargs...) confint(x; level=1-alpha, kwargs...)
\ No newline at end of file
diff --git a/src/durbin_watson.jl b/src/durbin_watson.jl
index 60e38fd7..3815f6d5 100644
--- a/src/durbin_watson.jl
+++ b/src/durbin_watson.jl
@@ -177,7 +177,7 @@ function pan_algorithm(a::AbstractArray, x::Float64, m::Int, n::Int)
 
 end
 
-function pvalue(x::DurbinWatsonTest; tail=:both)
+function StatsAPI.pvalue(x::DurbinWatsonTest; tail=:both)
 
     exact_problem_flag = 0
     if (x.p_compute == :ndep && x.n <= 100) || x.p_compute == :exact
diff --git a/src/f.jl b/src/f.jl
index d4f4e367..26dbae33 100644
--- a/src/f.jl
+++ b/src/f.jl
@@ -63,7 +63,7 @@ function show_params(io::IO, x::VarianceFTest, ident)
     println(io, ident, "degrees of freedom:     [$(x.df_x), $(x.df_y)]")
 end
 
-function pvalue(x::VarianceFTest; tail=:both)
+function StatsAPI.pvalue(x::VarianceFTest; tail=:both)
     dist = FDist(x.df_x, x.df_y)
     if tail == :both
         return 1 - 2*abs(cdf(dist, x.F) - 0.5)
diff --git a/src/fisher.jl b/src/fisher.jl
index 9b7a8ced..bde4e673 100644
--- a/src/fisher.jl
+++ b/src/fisher.jl
@@ -120,7 +120,7 @@ For `tail = :both`, possible values for `method` are:
     Blaker’s exact tests". Biostatistics, Volume 11, Issue 2, 1 April 2010, Pages 373–374,
     [link](https://doi.org/10.1093/biostatistics/kxp050)
 """
-function pvalue(x::FisherExactTest; tail=:both, method=:central)
+function StatsAPI.pvalue(x::FisherExactTest; tail=:both, method=:central)
     if tail == :both && method != :central
         if method == :minlike
             p = pvalue_both_minlike(x)
@@ -179,7 +179,7 @@ Fisher's non-central hypergeometric distribution. For `tail = :both`, the only
     Blaker’s exact tests". Biostatistics, Volume 11, Issue 2, 1 April 2010, Pages 373–374,
     [link](https://doi.org/10.1093/biostatistics/kxp050)
 """
-function StatsBase.confint(x::FisherExactTest; level::Float64=0.95, tail=:both, method=:central)
+function StatsAPI.confint(x::FisherExactTest; level::Float64=0.95, tail=:both, method=:central)
     check_level(level)
     if x.a == x.c == 0 || x.b == x.d == 0
         return (0.0, Inf)
@@ -203,8 +203,8 @@ function StatsBase.confint(x::FisherExactTest; level::Float64=0.95, tail=:both,
         end
     elseif tail == :both
         if method == :central
-            (StatsBase.confint(x, level=1-(1-level)/2, tail=:right)[1],
-             StatsBase.confint(x, level=1-(1-level)/2, tail=:left)[2])
+            (confint(x, level=1-(1-level)/2, tail=:right)[1],
+             confint(x, level=1-(1-level)/2, tail=:left)[2])
         else
             throw(ArgumentError("method=$(method) is not implemented yet"))
         end
diff --git a/src/hotelling.jl b/src/hotelling.jl
index 47635664..384eacac 100644
--- a/src/hotelling.jl
+++ b/src/hotelling.jl
@@ -5,7 +5,7 @@ export OneSampleHotellingT2Test, EqualCovHotellingT2Test, UnequalCovHotellingT2T
 abstract type HotellingT2TestTest <: HypothesisTest end
 
 default_tail(::HotellingT2TestTest) = :right
-pvalue(T::HotellingT2TestTest; tail=:right) = pvalue(FDist(dof(T)...), T.F, tail=tail)
+StatsAPI.pvalue(T::HotellingT2TestTest; tail=:right) = pvalue(FDist(dof(T)...), T.F, tail=tail)
 
 function show_params(io::IO, T::HotellingT2TestTest, indent="")
     println(io, indent, "number of observations: ", nobs(T))
@@ -42,8 +42,8 @@ struct OneSampleHotellingT2Test <: HotellingT2TestTest
     S::Matrix
 end
 
-StatsBase.nobs(T::OneSampleHotellingT2Test) = T.n
-StatsBase.dof(T::OneSampleHotellingT2Test) = (T.p, T.n - T.p)
+StatsAPI.nobs(T::OneSampleHotellingT2Test) = T.n
+StatsAPI.dof(T::OneSampleHotellingT2Test) = (T.p, T.n - T.p)
 
 """
     OneSampleHotellingT2Test(X::AbstractMatrix, μ₀=<zero vector>)
@@ -112,8 +112,8 @@ function EqualCovHotellingT2Test(X::AbstractMatrix, Y::AbstractMatrix)
     return EqualCovHotellingT2Test(T², F, nx, ny, p, Δ, S)
 end
 
-StatsBase.nobs(T::EqualCovHotellingT2Test) = (T.nx, T.ny)
-StatsBase.dof(T::EqualCovHotellingT2Test) = (T.p, T.nx + T.ny - T.p - 1)
+StatsAPI.nobs(T::EqualCovHotellingT2Test) = (T.nx, T.ny)
+StatsAPI.dof(T::EqualCovHotellingT2Test) = (T.p, T.nx + T.ny - T.p - 1)
 
 testname(::EqualCovHotellingT2Test) =
     "Two sample Hotelling's T² test (equal covariance matrices)"
@@ -156,8 +156,8 @@ function UnequalCovHotellingT2Test(X::AbstractMatrix, Y::AbstractMatrix)
     return UnequalCovHotellingT2Test(T², F, nx, ny, p, ν, Δ, ST)
 end
 
-StatsBase.nobs(T::UnequalCovHotellingT2Test) = (T.nx, T.ny)
-StatsBase.dof(T::UnequalCovHotellingT2Test) = (T.p, T.ν)
+StatsAPI.nobs(T::UnequalCovHotellingT2Test) = (T.nx, T.ny)
+StatsAPI.dof(T::UnequalCovHotellingT2Test) = (T.p, T.ν)
 
 testname(::UnequalCovHotellingT2Test) =
     "Two sample Hotelling's T² test (unequal covariance matrices)"
diff --git a/src/jarque_bera.jl b/src/jarque_bera.jl
index 15cb5f72..284798db 100644
--- a/src/jarque_bera.jl
+++ b/src/jarque_bera.jl
@@ -103,4 +103,4 @@ function show_params(io::IO, x::JarqueBeraTest, ident)
     println(io, ident, "JB statistic:                   ", x.JB)
 end
 
-pvalue(x::JarqueBeraTest) = pvalue(Chisq(2), x.JB; tail=:right)
+StatsAPI.pvalue(x::JarqueBeraTest) = pvalue(Chisq(2), x.JB; tail=:right)
diff --git a/src/kolmogorov_smirnov.jl b/src/kolmogorov_smirnov.jl
index 0a9dc776..2413c389 100644
--- a/src/kolmogorov_smirnov.jl
+++ b/src/kolmogorov_smirnov.jl
@@ -77,7 +77,7 @@ function show_params(io::IO, x::ExactOneSampleKSTest, ident="")
     println(io, ident, "number of observations:   $(x.n)")
 end
 
-function pvalue(x::ExactKSTest; tail=:both)
+function StatsAPI.pvalue(x::ExactKSTest; tail=:both)
     if tail == :left
         pvalue(KSOneSided(x.n), x.δn; tail=:right)
     elseif tail == :right
@@ -123,7 +123,7 @@ function show_params(io::IO, x::ApproximateOneSampleKSTest, ident="")
 end
 
 # one-sided: http://www.encyclopediaofmath.org/index.php/Kolmogorov-Smirnov_test
-function pvalue(x::ApproximateOneSampleKSTest; tail=:both)
+function StatsAPI.pvalue(x::ApproximateOneSampleKSTest; tail=:both)
     if tail == :left
         exp(-2*x.n*x.δn^2)
     elseif tail == :right
@@ -175,7 +175,7 @@ function show_params(io::IO, x::ApproximateTwoSampleKSTest, ident="")
     println(io, ident, "KS-statistic:              $(sqrt(n)*x.δ)")
 end
 
-function pvalue(x::ApproximateTwoSampleKSTest; tail=:both)
+function StatsAPI.pvalue(x::ApproximateTwoSampleKSTest; tail=:both)
     n = x.n_x*x.n_y/(x.n_x+x.n_y)
     if tail == :left
         exp(-2*n*x.δn^2)
diff --git a/src/kruskal_wallis.jl b/src/kruskal_wallis.jl
index 48e80a19..934b6866 100644
--- a/src/kruskal_wallis.jl
+++ b/src/kruskal_wallis.jl
@@ -92,7 +92,7 @@ function show_params(io::IO, x::KruskalWallisTest, ident)
     println(io, ident, "adjustment for ties:                 ", x.tie_adjustment)
 end
 
-pvalue(x::KruskalWallisTest) = pvalue(Chisq(x.df), x.H; tail=:right)
+StatsAPI.pvalue(x::KruskalWallisTest) = pvalue(Chisq(x.df), x.H; tail=:right)
 
 
 ## helper
diff --git a/src/mann_whitney.jl b/src/mann_whitney.jl
index 1bcd38b3..2422dacb 100644
--- a/src/mann_whitney.jl
+++ b/src/mann_whitney.jl
@@ -133,7 +133,7 @@ function mwuenumerate(x::ExactMannWhitneyUTest)
     (le/tot, gr/tot)
 end
 
-function pvalue(x::ExactMannWhitneyUTest; tail=:both)
+function StatsAPI.pvalue(x::ExactMannWhitneyUTest; tail=:both)
     check_tail(tail)
 
     if x.tie_adjustment == 0
@@ -221,7 +221,7 @@ function show_params(io::IO, x::ApproximateMannWhitneyUTest, ident)
     println(io, ident, "normal approximation (μ, σ):          ", (x.mu, x.sigma))
 end
 
-function pvalue(x::ApproximateMannWhitneyUTest; tail=:both)
+function StatsAPI.pvalue(x::ApproximateMannWhitneyUTest; tail=:both)
     check_tail(tail)
 
     if x.mu == x.sigma == 0
diff --git a/src/permutation.jl b/src/permutation.jl
index 918dfb7b..5242db21 100644
--- a/src/permutation.jl
+++ b/src/permutation.jl
@@ -44,7 +44,7 @@ ApproximatePermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
                            f::Function, n::Int) where {R<:Real,S<:Real} =
     ApproximatePermutationTest(Random.default_rng(), x, y, f, n)
 
-function pvalue(apt::PermutationTest; tail=:both)
+function StatsAPI.pvalue(apt::PermutationTest; tail=:both)
     if tail == :both
         count = sum(abs(apt.observation) <= abs(x) for x in apt.samples)
     elseif tail == :left
diff --git a/src/power_divergence.jl b/src/power_divergence.jl
index 39f6b844..5cfa7949 100644
--- a/src/power_divergence.jl
+++ b/src/power_divergence.jl
@@ -41,7 +41,7 @@ end
 population_param_of_interest(x::PowerDivergenceTest) = ("Multinomial Probabilities", x.theta0, x.thetahat)
 default_tail(test::PowerDivergenceTest) = :right
 
-pvalue(x::PowerDivergenceTest; tail=:right) = pvalue(Chisq(x.df),x.stat; tail=tail)
+StatsAPI.pvalue(x::PowerDivergenceTest; tail=:right) = pvalue(Chisq(x.df),x.stat; tail=tail)
 
 """
     confint(test::PowerDivergenceTest; level = 0.95, tail = :both, method = :auto)
@@ -67,18 +67,18 @@ one of the following methods. Possible values for `method` are:
   * Gold, R. Z. Tests Auxiliary to ``χ^2`` Tests in a Markov Chain. Annals of
     Mathematical Statistics, 30:56-74, 1963.
 """
-function StatsBase.confint(x::PowerDivergenceTest; level::Float64=0.95,
-                           tail::Symbol=:both, method::Symbol=:auto, correct::Bool=true,
-                           bootstrap_iters::Int64=10000, GC::Bool=true)
+function StatsAPI.confint(x::PowerDivergenceTest; level::Float64=0.95,
+                          tail::Symbol=:both, method::Symbol=:auto, correct::Bool=true,
+                          bootstrap_iters::Int64=10000, GC::Bool=true)
     check_level(level)
 
     m  = length(x.thetahat)
 
     if tail == :left
-        i = StatsBase.confint(x, level=1-(1-level)*2, method=method, GC=GC)
+        i = confint(x, level=1-(1-level)*2, method=method, GC=GC)
         Tuple{Float64,Float64}[(0.0, i[j][2]) for j in 1:m]
     elseif tail == :right
-        i = StatsBase.confint(x, level=1-(1-level)*2, method=method, GC=GC)
+        i = confint(x, level=1-(1-level)*2, method=method, GC=GC)
         Tuple{Float64,Float64}[(i[j][1], 1.0) for j in 1:m]
     elseif tail == :both
         if method == :auto
diff --git a/src/t.jl b/src/t.jl
index d2896aa4..0ed1c6df 100644
--- a/src/t.jl
+++ b/src/t.jl
@@ -28,18 +28,18 @@ export OneSampleTTest, TwoSampleTTest, EqualVarianceTTest,
 abstract type TTest <: HypothesisTest end
 abstract type TwoSampleTTest <: TTest end
 
-pvalue(x::TTest; tail=:both) = pvalue(TDist(x.df), x.t; tail=tail)
+StatsAPI.pvalue(x::TTest; tail=:both) = pvalue(TDist(x.df), x.t; tail=tail)
 
 default_tail(test::TTest) = :both
 
 # confidence interval by inversion
-function StatsBase.confint(x::TTest; level::Float64=0.95, tail=:both)
+function StatsAPI.confint(x::TTest; level::Float64=0.95, tail=:both)
     check_level(level)
 
     if tail == :left
-        (-Inf, StatsBase.confint(x, level=1-(1-level)*2)[2])
+        (-Inf, confint(x, level=1-(1-level)*2)[2])
     elseif tail == :right
-        (StatsBase.confint(x, level=1-(1-level)*2)[1], Inf)
+        (confint(x, level=1-(1-level)*2)[1], Inf)
     elseif tail == :both
         q = quantile(TDist(x.df), 1-(1-level)/2)
         (x.xbar-q*x.stderr, x.xbar+q*x.stderr)
diff --git a/src/var_equality.jl b/src/var_equality.jl
index c80a8d69..3e1775b1 100644
--- a/src/var_equality.jl
+++ b/src/var_equality.jl
@@ -18,16 +18,16 @@ testname(t::VarianceEqualityTest) = t.description[1]
 teststatisticname(t::VarianceEqualityTest{TD}) where {TD <: ContinuousDistribution} =
     length(t.description[3]) != 0 ? t.description[3] : (TD <: FDist ? "F" : "χ²")
 
-StatsBase.nobs(t::VarianceEqualityTest) = t.Nᵢ
-StatsBase.dof(t::VarianceEqualityTest{Chisq}) = t.DFt
-StatsBase.dof(t::VarianceEqualityTest{FDist}) = (t.DFt, t.DFe)
+StatsAPI.nobs(t::VarianceEqualityTest) = t.Nᵢ
+StatsAPI.dof(t::VarianceEqualityTest{Chisq}) = t.DFt
+StatsAPI.dof(t::VarianceEqualityTest{FDist}) = (t.DFt, t.DFe)
 
 teststatistic(t::VarianceEqualityTest{FDist}) = (t.DFe/t.DFt)*sum(t.SStᵢ)/sum(t.SSeᵢ)
 function teststatistic(t::VarianceEqualityTest{Chisq})
     y = sum(t.SStᵢ)/sum(t.SSeᵢ)
     y*(t.DFe+t.DFt)/(1 + y) # sum(t.SStᵢ)/t.s²
 end
-pvalue(t::VarianceEqualityTest{TD}; tail=:right) where {TD <: ContinuousDistribution} =
+StatsAPI.pvalue(t::VarianceEqualityTest{TD}; tail=:right) where {TD <: ContinuousDistribution} =
     pvalue(TD(dof(t)...), teststatistic(t), tail=tail)
 
 function show_params(io::IO, t::VarianceEqualityTest{TD},
@@ -49,11 +49,11 @@ function Base.show(io::IOContext, t::VarianceEqualityTest)
         MSe = SSe/t.DFe
         println(io, "Source            SS    DF        MS         F  P-value")
         println(io, repeat("-", 55))
-        StatsBase.@printf(io, "Treatments  %8.3f  %4d  %8.3f  %8.3f  %7.5f\n",
-                              SSt, t.DFt, MSt, MSt/MSe, pvalue(t))
-        StatsBase.@printf(io, "Error       %8.3f  %4d  %8.3f\n", SSe, t.DFe, MSe)
+        @printf(io, "Treatments  %8.3f  %4d  %8.3f  %8.3f  %7.5f\n",
+                SSt, t.DFt, MSt, MSt/MSe, pvalue(t))
+        @printf(io, "Error       %8.3f  %4d  %8.3f\n", SSe, t.DFe, MSe)
         println(io, repeat("-", 55))
-        StatsBase.@printf(io, "Total       %8.3f  %4d\n", SSt+SSe, t.DFt+t.DFe)
+        @printf(io, "Total       %8.3f  %4d\n", SSt+SSe, t.DFt+t.DFe)
     end
 end
 
diff --git a/src/wald_wolfowitz.jl b/src/wald_wolfowitz.jl
index c847f462..1843fe02 100644
--- a/src/wald_wolfowitz.jl
+++ b/src/wald_wolfowitz.jl
@@ -12,7 +12,7 @@ end
 testname(::WaldWolfowitzTest) = "Wald-Wolfowitz Test"
 population_param_of_interest(x::WaldWolfowitzTest) = ("Number of runs", x.μ, x.nruns) # parameter of interest: name, value under h0, point estimate
 default_tail(::WaldWolfowitzTest) = :both
-pvalue(test::WaldWolfowitzTest; tail=:both) = pvalue(Normal(0.0, 1.0), test.z; tail=tail)
+StatsAPI.pvalue(test::WaldWolfowitzTest; tail=:both) = pvalue(Normal(0.0, 1.0), test.z; tail=tail)
 
 
 function show_params(io::IO, x::WaldWolfowitzTest, ident="")
diff --git a/src/white.jl b/src/white.jl
index 93ba8a96..6af008a7 100644
--- a/src/white.jl
+++ b/src/white.jl
@@ -116,5 +116,5 @@ function show_params(io::IO, t::WhiteTest, ident = "")
     println(io, ident, "type:                   ", t.type)
 end
 
-StatsBase.dof(t::WhiteTest) = t.dof
-pvalue(t::WhiteTest) = pvalue(Chisq(t.dof), t.lm, tail=:right)
+StatsAPI.dof(t::WhiteTest) = t.dof
+StatsAPI.pvalue(t::WhiteTest) = pvalue(Chisq(t.dof), t.lm, tail=:right)
diff --git a/src/wilcoxon.jl b/src/wilcoxon.jl
index 5756165f..f3e5f793 100644
--- a/src/wilcoxon.jl
+++ b/src/wilcoxon.jl
@@ -132,7 +132,7 @@ function signedrankenumerate(x::ExactSignedRankTest)
     (le/tot, gr/tot)
 end
 
-function pvalue(x::ExactSignedRankTest; tail=:both)
+function StatsAPI.pvalue(x::ExactSignedRankTest; tail=:both)
     check_tail(tail)
 
     n = length(x.ranks)
@@ -163,7 +163,7 @@ function pvalue(x::ExactSignedRankTest; tail=:both)
     end
 end
 
-StatsBase.confint(x::ExactSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)
+StatsAPI.confint(x::ExactSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)
 
 
 ## APPROXIMATE SIGNED RANK TEST
@@ -224,7 +224,7 @@ function show_params(io::IO, x::ApproximateSignedRankTest, ident)
     println(io, ident, "normal approximation (μ, σ): ", (x.mu, x.sigma))
 end
 
-function pvalue(x::ApproximateSignedRankTest; tail=:both)
+function StatsAPI.pvalue(x::ApproximateSignedRankTest; tail=:both)
     check_tail(tail)
 
     if x.mu == x.sigma == 0
@@ -238,7 +238,7 @@ function pvalue(x::ApproximateSignedRankTest; tail=:both)
     end
 end
 
-StatsBase.confint(x::ApproximateSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)
+StatsAPI.confint(x::ApproximateSignedRankTest; level::Real=0.95, tail=:both) = calculate_ci(x.vals, level, tail=tail)
 
 # implementation method inspired by these notes: http://www.stat.umn.edu/geyer/old03/5102/notes/rank.pdf
 function calculate_ci(x::AbstractVector, level::Real=0.95; tail=:both)
diff --git a/src/z.jl b/src/z.jl
index af3370ef..141d7da2 100644
--- a/src/z.jl
+++ b/src/z.jl
@@ -28,18 +28,18 @@ export OneSampleZTest, TwoSampleZTest, EqualVarianceZTest,
 abstract type ZTest <: HypothesisTest end
 abstract type TwoSampleZTest <: ZTest end
 
-pvalue(x::ZTest; tail=:both) = pvalue(Normal(0.0, 1.0), x.z; tail=tail)
+StatsAPI.pvalue(x::ZTest; tail=:both) = pvalue(Normal(0.0, 1.0), x.z; tail=tail)
 
 default_tail(test::ZTest) = :both
 
 # confidence interval by inversion
-function StatsBase.confint(x::ZTest; level::Float64=0.95, tail=:both)
+function StatsAPI.confint(x::ZTest; level::Float64=0.95, tail=:both)
     check_level(level)
 
     if tail == :left
-        (-Inf, StatsBase.confint(x, level=1-(1-level)*2)[2])
+        (-Inf, confint(x, level=1-(1-level)*2)[2])
     elseif tail == :right
-        (StatsBase.confint(x, level=1-(1-level)*2)[1], Inf)
+        (confint(x, level=1-(1-level)*2)[1], Inf)
     elseif tail == :both
         q = cquantile(Normal(0.0, 1.0), (1-level)/2)
         (x.xbar-q*x.stderr, x.xbar+q*x.stderr)
diff --git a/test/bartlett.jl b/test/bartlett.jl
index 9eaefd89..612562ec 100644
--- a/test/bartlett.jl
+++ b/test/bartlett.jl
@@ -1,6 +1,5 @@
 using HypothesisTests
 using Test
-using StatsBase
 using DelimitedFiles
 
 @testset "Bartlett's test" begin
diff --git a/test/correlation.jl b/test/correlation.jl
index d2686378..a8c10c3e 100644
--- a/test/correlation.jl
+++ b/test/correlation.jl
@@ -1,7 +1,6 @@
 using HypothesisTests
 using Test
 using DelimitedFiles
-using StatsBase
 
 @testset "Correlation" begin
     # Columns are line number, calcium, iron
diff --git a/test/var_equality.jl b/test/var_equality.jl
index e58383ce..604ba947 100644
--- a/test/var_equality.jl
+++ b/test/var_equality.jl
@@ -1,6 +1,5 @@
 using HypothesisTests
 using Test
-using StatsBase
 using Statistics
 using DelimitedFiles