diff --git a/REQUIRE b/REQUIRE
index b1ef85825..ae6850ade 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,4 +1,4 @@
-julia 0.7.0-alpha
+julia 0.7.0-beta2
 DataStructures 0.5.0
 SortingAlgorithms
 Missings
diff --git a/appveyor.yml b/appveyor.yml
index 9b6671c60..cef728064 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -1,9 +1,17 @@
 environment:
   matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.7/julia-0.7-latest-win32.exe"
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.7/julia-0.7-latest-win64.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
-  - JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
+  - julia_version: 0.7
+  - julia_version: latest
+
+platform:
+  - x86
+  - x64
+
+## uncomment the following lines to allow failures on nightly julia
+## (tests will run but not make your overall status red)
+#matrix:
+#  allow_failures:
+#  - julia_version: latest
 
 branches:
   only:
@@ -17,24 +25,12 @@ notifications:
     on_build_status_changed: false
 
 install:
-  - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
-# If there's a newer build queued for the same PR, cancel this one
-  - ps: if ($env:APPVEYOR_PULL_REQUEST_NUMBER -and $env:APPVEYOR_BUILD_NUMBER -ne ((Invoke-RestMethod `
-        https://ci.appveyor.com/api/projects/$env:APPVEYOR_ACCOUNT_NAME/$env:APPVEYOR_PROJECT_SLUG/history?recordsNumber=50).builds | `
-        Where-Object pullRequestId -eq $env:APPVEYOR_PULL_REQUEST_NUMBER)[0].buildNumber) { `
-        throw "There are newer queued builds for this pull request, failing early." }
-# Download most recent Julia Windows binary
-  - ps: (new-object net.webclient).DownloadFile(
-        $env:JULIA_URL,
-        "C:\projects\julia-binary.exe")
-# Run installer silently, output to C:\projects\julia
-  - C:\projects\julia-binary.exe /S /D=C:\projects\julia
+  - ps: iex ((new-object net.webclient).DownloadString('https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/master/bin/install.ps1'))
 
 build_script:
-# Need to convert from shallow to complete for Pkg.clone to work
-  - IF EXIST .git\shallow (git fetch --unshallow)
-  - C:\projects\julia\bin\julia -e "versioninfo();
-      Pkg.clone(pwd(), \"StatsBase\"); Pkg.build(\"StatsBase\")"
+  - echo "%JL_BUILD_SCRIPT%"
+  - julia -e "%JL_BUILD_SCRIPT%"
 
 test_script:
-  - C:\projects\julia\bin\julia -e "Pkg.test(\"StatsBase\")"
+  - echo "%JL_TEST_SCRIPT%"
+  - julia -e "%JL_TEST_SCRIPT%"
diff --git a/src/StatsBase.jl b/src/StatsBase.jl
index 8926787b8..8b5873396 100644
--- a/src/StatsBase.jl
+++ b/src/StatsBase.jl
@@ -8,27 +8,16 @@ module StatsBase
     using SortingAlgorithms
     using Missings
 
-    if VERSION >= v"0.7.0-beta.85"
-        using Statistics
-        import Statistics: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
-                           cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
-                           median, middle
-    else
-        import Base: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
-                     cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
-                     median, middle
-    end
-
-    # Location of the `mean` function, used for calling the function inside of
-    # functions definitions that have `mean` as a keyword argument
-    const MEANHOME = (VERSION >= v"0.7.0-beta.85" ? Statistics : Base)::Module
-
+    using Statistics
     using LinearAlgebra
     using Random
     using Printf
     using SparseArrays
     import Random: rand, rand!
     import LinearAlgebra: BlasReal, BlasFloat
+    import Statistics: mean, mean!, var, varm, varm!, std, stdm, cov, covm,
+                       cor, corm, cov2cor!, unscaled_covzm, quantile, sqrt!,
+                       median, middle
 
     ## tackle compatibility issues
 
@@ -206,18 +195,6 @@ module StatsBase
 
     ConvergenceException
 
-const BASESTATS_IN_STATSBASE = v"0.7.0-DEV.5238" <= VERSION < v"0.7.0-beta.85"
-
-@static if BASESTATS_IN_STATSBASE
-    export cor, cov, std, stdm, var, varm, linreg
-    include("base.jl")
-end
-
-module StatsCompat
-    import ..StatsBase: var, std, varm, cov, cor, MEANHOME
-end
-export StatsCompat
-
     # source files
 
     include("common.jl")
diff --git a/src/base.jl b/src/base.jl
deleted file mode 100644
index eea479415..000000000
--- a/src/base.jl
+++ /dev/null
@@ -1,617 +0,0 @@
-# This file contains code formerly part of Julia. License is MIT: https://julialang.org/license
-
-# deprecations from base/deprecated.jl
-@deprecate varm(A::AbstractArray, m::AbstractArray, dims; kwargs...) varm(A, m; kwargs..., dims=dims)
-@deprecate var(A::AbstractArray, dims; kwargs...)                    var(A; kwargs..., dims=dims)
-@deprecate std(A::AbstractArray, dims; kwargs...)                    std(A; kwargs..., dims=dims)
-@deprecate cov(X::AbstractMatrix, dim::Int; kwargs...)               cov(X; kwargs..., dims=dim)
-@deprecate cov(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int; kwargs...) cov(x, y; kwargs..., dims=dim)
-@deprecate cor(X::AbstractMatrix, dim::Int)                          cor(X, dims=dim)
-@deprecate cor(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int)   cor(x, y, dims=dim)
-
-# PR #21709
-@deprecate cov(x::AbstractVector, corrected::Bool) cov(x, corrected=corrected)
-@deprecate cov(x::AbstractMatrix, vardim::Int, corrected::Bool) cov(x, dims=vardim, corrected=corrected)
-@deprecate cov(X::AbstractVector, Y::AbstractVector, corrected::Bool) cov(X, Y, corrected=corrected)
-@deprecate cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int, corrected::Bool) cov(X, Y, dims=vardim, corrected=corrected)
-
-##### variances #####
-
-# faster computation of real(conj(x)*y)
-realXcY(x::Real, y::Real) = x*y
-realXcY(x::Complex, y::Complex) = real(x)*real(y) + imag(x)*imag(y)
-
-var(iterable; corrected::Bool=true, mean=nothing) = _var(iterable, corrected, mean)
-
-function _var(iterable, corrected::Bool, mean)
-    y = iterate(iterable)
-    if y === nothing
-        throw(ArgumentError("variance of empty collection undefined: $(repr(iterable))"))
-    end
-    count = 1
-    value, state = y
-    y = iterate(iterable, state)
-    if mean === nothing
-        # Use Welford algorithm as seen in (among other places)
-        # Knuth's TAOCP, Vol 2, page 232, 3rd edition.
-        M = value / 1
-        S = real(zero(M))
-        while y !== nothing
-            value, state = y
-            y = iterate(iterable, state)
-            count += 1
-            new_M = M + (value - M) / count
-            S = S + realXcY(value - M, value - new_M)
-            M = new_M
-        end
-        return S / (count - Int(corrected))
-    elseif isa(mean, Number) # mean provided
-        # Cannot use a compensated version, e.g. the one from
-        # "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances."
-        # by Chan, Golub, and LeVeque, Technical Report STAN-CS-79-773,
-        # Department of Computer Science, Stanford University,
-        # because user can provide mean value that is different to mean(iterable)
-        sum2 = abs2(value - mean::Number)
-        while y !== nothing
-            value, state = y
-            y = iterate(iterable, state)
-            count += 1
-            sum2 += abs2(value - mean)
-        end
-        return sum2 / (count - Int(corrected))
-    else
-        throw(ArgumentError("invalid value of mean, $(mean)::$(typeof(mean))"))
-    end
-end
-
-centralizedabs2fun(m) = x -> abs2.(x - m)
-centralize_sumabs2(A::AbstractArray, m) =
-    mapreduce(centralizedabs2fun(m), +, A)
-centralize_sumabs2(A::AbstractArray, m, ifirst::Int, ilast::Int) =
-    Base.mapreduce_impl(centralizedabs2fun(m), +, A, ifirst, ilast)
-
-function centralize_sumabs2!(R::AbstractArray{S}, A::AbstractArray, means::AbstractArray) where S
-    # following the implementation of _mapreducedim! at base/reducedim.jl
-    lsiz = Base.check_reducedims(R,A)
-    isempty(R) || fill!(R, zero(S))
-    isempty(A) && return R
-
-    if Base.has_fast_linear_indexing(A) && lsiz > 16
-        nslices = div(Base._length(A), lsiz)
-        ibase = first(LinearIndices(A))-1
-        for i = 1:nslices
-            @inbounds R[i] = centralize_sumabs2(A, means[i], ibase+1, ibase+lsiz)
-            ibase += lsiz
-        end
-        return R
-    end
-    indsAt, indsRt = Base.safe_tail(axes(A)), Base.safe_tail(axes(R)) # handle d=1 manually
-    keep, Idefault = Broadcast.shapeindexer(indsRt)
-    if Base.reducedim1(R, A)
-        i1 = first(Base.indices1(R))
-        @inbounds for IA in CartesianIndices(indsAt)
-            IR = Broadcast.newindex(IA, keep, Idefault)
-            r = R[i1,IR]
-            m = means[i1,IR]
-            @simd for i in axes(A, 1)
-                r += abs2(A[i,IA] - m)
-            end
-            R[i1,IR] = r
-        end
-    else
-        @inbounds for IA in CartesianIndices(indsAt)
-            IR = Broadcast.newindex(IA, keep, Idefault)
-            @simd for i in axes(A, 1)
-                R[i,IR] += abs2(A[i,IA] - means[i,IR])
-            end
-        end
-    end
-    return R
-end
-
-function varm!(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; corrected::Bool=true) where S
-    if isempty(A)
-        fill!(R, convert(S, NaN))
-    else
-        rn = div(Base._length(A), Base._length(R)) - Int(corrected)
-        centralize_sumabs2!(R, A, m)
-        R .= R .* (1 // rn)
-    end
-    return R
-end
-
-"""
-    varm(v, m; dims, corrected::Bool=true)
-
-Compute the sample variance of a collection `v` with known mean(s) `m`,
-optionally over the given dimensions. `m` may contain means for each dimension of
-`v`. If `corrected` is `true`, then the sum is scaled with `n-1`,
-whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`.
-
-!!! note
-    Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
-    to represent missing values, and the [`skipmissing`](@ref) function to omit them.
-"""
-varm(A::AbstractArray, m::AbstractArray; corrected::Bool=true, dims=:) = _varm(A, m, corrected, dims)
-
-_varm(A::AbstractArray{T}, m, corrected::Bool, region) where {T} =
-    varm!(Base.reducedim_init(t -> abs2(t)/2, +, A, region), A, m; corrected=corrected)
-
-varm(A::AbstractArray, m; corrected::Bool=true) = _varm(A, m, corrected, :)
-
-function _varm(A::AbstractArray{T}, m, corrected::Bool, ::Colon) where T
-    n = Base._length(A)
-    n == 0 && return typeof((abs2(zero(T)) + abs2(zero(T)))/2)(NaN)
-    return centralize_sumabs2(A, m) / (n - Int(corrected))
-end
-
-
-"""
-    var(v; dims, corrected::Bool=true, mean=nothing)
-
-Compute the sample variance of a vector or array `v`, optionally along the given dimensions.
-The algorithm will return an estimator of the generative distribution's variance
-under the assumption that each entry of `v` is an IID drawn from that generative
-distribution. This computation is equivalent to calculating `sum(abs2, v - mean(v)) /
-(length(v) - 1)`. If `corrected` is `true`, then the sum is scaled with `n-1`,
-whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`.
-The mean `mean` over the region may be provided.
-
-!!! note
-    Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
-    to represent missing values, and the [`skipmissing`](@ref) function to omit them.
-"""
-var(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _var(A, corrected, mean, dims)
-
-_var(A::AbstractArray, corrected::Bool, mean, dims) =
-    varm(A, something(mean, MEANHOME.mean(A, dims=dims)); corrected=corrected, dims=dims)
-
-_var(A::AbstractArray, corrected::Bool, mean, ::Colon) =
-    real(varm(A, something(mean, MEANHOME.mean(A)); corrected=corrected))
-
-varm(iterable, m; corrected::Bool=true) = _var(iterable, corrected, m)
-
-## variances over ranges
-
-varm(v::AbstractRange, m::AbstractArray) = range_varm(v, m)
-varm(v::AbstractRange, m) = range_varm(v, m)
-
-function range_varm(v::AbstractRange, m)
-    f  = first(v) - m
-    s  = step(v)
-    l  = length(v)
-    vv = f^2 * l / (l - 1) + f * s * l + s^2 * l * (2 * l - 1) / 6
-    if l == 0 || l == 1
-        return typeof(vv)(NaN)
-    end
-    return vv
-end
-
-function var(v::AbstractRange)
-    s  = step(v)
-    l  = length(v)
-    vv = abs2(s) * (l + 1) * l / 12
-    if l == 0 || l == 1
-        return typeof(vv)(NaN)
-    end
-    return vv
-end
-
-
-##### standard deviation #####
-
-function sqrt!(A::AbstractArray)
-    for i in eachindex(A)
-        @inbounds A[i] = sqrt(A[i])
-    end
-    A
-end
-
-stdm(A::AbstractArray, m; corrected::Bool=true) =
-    sqrt.(varm(A, m; corrected=corrected))
-
-"""
-    std(v; corrected::Bool=true, mean=nothing, dims)
-
-Compute the sample standard deviation of a vector or array `v`, optionally along the given
-dimensions. The algorithm returns an estimator of the generative distribution's standard
-deviation under the assumption that each entry of `v` is an IID drawn from that generative
-distribution. This computation is equivalent to calculating `sqrt(sum((v - mean(v)).^2) /
-(length(v) - 1))`. A pre-computed `mean` may be provided. If `corrected` is `true`,
-then the sum is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is
-`false` where `n = length(x)`.
-
-!!! note
-    Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
-    to represent missing values, and the [`skipmissing`](@ref) function to omit them.
-"""
-std(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _std(A, corrected, mean, dims)
-
-_std(A::AbstractArray, corrected::Bool, mean, dims) =
-    sqrt.(var(A; corrected=corrected, mean=mean, dims=dims))
-
-_std(A::AbstractArray, corrected::Bool, mean, ::Colon) =
-    sqrt.(var(A; corrected=corrected, mean=mean))
-
-_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, dims) =
-    sqrt!(var(A; corrected=corrected, mean=mean, dims=dims))
-
-_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, ::Colon) =
-    sqrt.(var(A; corrected=corrected, mean=mean))
-
-std(iterable; corrected::Bool=true, mean=nothing) =
-    sqrt(var(iterable, corrected=corrected, mean=mean))
-
-"""
-    stdm(v, m; corrected::Bool=true)
-
-Compute the sample standard deviation of a vector `v`
-with known mean `m`. If `corrected` is `true`,
-then the sum is scaled with `n-1`, whereas the sum is
-scaled with `n` if `corrected` is `false` where `n = length(x)`.
-
-!!! note
-    Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type
-    to represent missing values, and the [`skipmissing`](@ref) function to omit them.
-"""
-stdm(iterable, m; corrected::Bool=true) =
-    std(iterable, corrected=corrected, mean=m)
-
-
-###### covariance ######
-
-# auxiliary functions
-
-_conj(x::AbstractArray{<:Real}) = x
-_conj(x::AbstractArray) = conj(x)
-
-_getnobs(x::AbstractVector, vardim::Int) = Base._length(x)
-_getnobs(x::AbstractMatrix, vardim::Int) = size(x, vardim)
-
-function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int)
-    n = _getnobs(x, vardim)
-    _getnobs(y, vardim) == n || throw(DimensionMismatch("dimensions of x and y mismatch"))
-    return n
-end
-
-_vmean(x::AbstractVector, vardim::Int) = mean(x)
-_vmean(x::AbstractMatrix, vardim::Int) = mean(x, dims=vardim)
-
-# core functions
-
-unscaled_covzm(x::AbstractVector{<:Number})    = sum(abs2, x)
-unscaled_covzm(x::AbstractVector)              = sum(t -> t*t', x)
-unscaled_covzm(x::AbstractMatrix, vardim::Int) = (vardim == 1 ? _conj(x'x) : x * x')
-
-unscaled_covzm(x::AbstractVector, y::AbstractVector) = dot(y, x)
-unscaled_covzm(x::AbstractVector, y::AbstractMatrix, vardim::Int) =
-    (vardim == 1 ? *(transpose(x), _conj(y)) : *(transpose(x), transpose(_conj(y))))
-unscaled_covzm(x::AbstractMatrix, y::AbstractVector, vardim::Int) =
-    (c = vardim == 1 ? *(transpose(x), _conj(y)) :  x * _conj(y); reshape(c, length(c), 1))
-unscaled_covzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int) =
-    (vardim == 1 ? *(transpose(x), _conj(y)) : *(x, adjoint(y)))
-
-# covzm (with centered data)
-
-covzm(x::AbstractVector; corrected::Bool=true) = unscaled_covzm(x) / (Base._length(x) - Int(corrected))
-function covzm(x::AbstractMatrix, vardim::Int=1; corrected::Bool=true)
-    C = unscaled_covzm(x, vardim)
-    T = promote_type(typeof(first(C) / 1), eltype(C))
-    A = convert(AbstractMatrix{T}, C)
-    b = 1//(size(x, vardim) - corrected)
-    A .= A .* b
-    return A
-end
-covzm(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
-    unscaled_covzm(x, y) / (Base._length(x) - Int(corrected))
-function covzm(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int=1; corrected::Bool=true)
-    C = unscaled_covzm(x, y, vardim)
-    T = promote_type(typeof(first(C) / 1), eltype(C))
-    A = convert(AbstractArray{T}, C)
-    b = 1//(_getnobs(x, y, vardim) - corrected)
-    A .= A .* b
-    return A
-end
-
-# covm (with provided mean)
-## Use map(t -> t - xmean, x) instead of x .- xmean to allow for Vector{Vector}
-## which can't be handled by broadcast
-covm(x::AbstractVector, xmean; corrected::Bool=true) =
-    covzm(map(t -> t - xmean, x); corrected=corrected)
-covm(x::AbstractMatrix, xmean, vardim::Int=1; corrected::Bool=true) =
-    covzm(x .- xmean, vardim; corrected=corrected)
-covm(x::AbstractVector, xmean, y::AbstractVector, ymean; corrected::Bool=true) =
-    covzm(map(t -> t - xmean, x), map(t -> t - ymean, y); corrected=corrected)
-covm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean, vardim::Int=1; corrected::Bool=true) =
-    covzm(x .- xmean, y .- ymean, vardim; corrected=corrected)
-
-# cov (API)
-"""
-    cov(x::AbstractVector; corrected::Bool=true)
-
-Compute the variance of the vector `x`. If `corrected` is `true` (the default) then the sum
-is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`.
-"""
-cov(x::AbstractVector; corrected::Bool=true) = covm(x, mean(x); corrected=corrected)
-
-"""
-    cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
-
-Compute the covariance matrix of the matrix `X` along the dimension `dims`. If `corrected`
-is `true` (the default) then the sum is scaled with `n-1`, whereas the sum is scaled with `n`
-if `corrected` is `false` where `n = size(X, dims)`.
-"""
-cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true) =
-    covm(X, _vmean(X, dims), dims; corrected=corrected)
-
-"""
-    cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
-
-Compute the covariance between the vectors `x` and `y`. If `corrected` is `true` (the
-default), computes ``\\frac{1}{n-1}\\sum_{i=1}^n (x_i-\\bar x) (y_i-\\bar y)^*`` where
-``*`` denotes the complex conjugate and `n = length(x) = length(y)`. If `corrected` is
-`false`, computes ``\\frac{1}{n}\\sum_{i=1}^n (x_i-\\bar x) (y_i-\\bar y)^*``.
-"""
-cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
-    covm(x, mean(x), y, mean(y); corrected=corrected)
-
-"""
-    cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true)
-
-Compute the covariance between the vectors or matrices `X` and `Y` along the dimension
-`dims`. If `corrected` is `true` (the default) then the sum is scaled with `n-1`, whereas
-the sum is scaled with `n` if `corrected` is `false` where `n = size(X, dims) = size(Y, dims)`.
-"""
-cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true) =
-    covm(X, _vmean(X, dims), Y, _vmean(Y, dims), dims; corrected=corrected)
-
-##### correlation #####
-
-"""
-    clampcor(x)
-
-Clamp a real correlation to between -1 and 1, leaving complex correlations unchanged
-"""
-clampcor(x::Real) = clamp(x, -1, 1)
-clampcor(x) = x
-
-# cov2cor!
-
-function cov2cor!(C::AbstractMatrix{T}, xsd::AbstractArray) where T
-    nx = length(xsd)
-    size(C) == (nx, nx) || throw(DimensionMismatch("inconsistent dimensions"))
-    for j = 1:nx
-        for i = 1:j-1
-            C[i,j] = adjoint(C[j,i])
-        end
-        C[j,j] = oneunit(T)
-        for i = j+1:nx
-            C[i,j] = clampcor(C[i,j] / (xsd[i] * xsd[j]))
-        end
-    end
-    return C
-end
-function cov2cor!(C::AbstractMatrix, xsd, ysd::AbstractArray)
-    nx, ny = size(C)
-    length(ysd) == ny || throw(DimensionMismatch("inconsistent dimensions"))
-    for (j, y) in enumerate(ysd)   # fixme (iter): here and in all `cov2cor!` we assume that `C` is efficiently indexed by integers
-        for i in 1:nx
-            C[i,j] = clampcor(C[i, j] / (xsd * y))
-        end
-    end
-    return C
-end
-function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd)
-    nx, ny = size(C)
-    length(xsd) == nx || throw(DimensionMismatch("inconsistent dimensions"))
-    for j in 1:ny
-        for (i, x) in enumerate(xsd)
-            C[i,j] = clampcor(C[i,j] / (x * ysd))
-        end
-    end
-    return C
-end
-function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray)
-    nx, ny = size(C)
-    (length(xsd) == nx && length(ysd) == ny) ||
-        throw(DimensionMismatch("inconsistent dimensions"))
-    for (i, x) in enumerate(xsd)
-        for (j, y) in enumerate(ysd)
-            C[i,j] = clampcor(C[i,j] / (x * y))
-        end
-    end
-    return C
-end
-
-# corzm (non-exported, with centered data)
-
-corzm(x::AbstractVector{T}) where {T} = one(real(T))
-function corzm(x::AbstractMatrix, vardim::Int=1)
-    c = unscaled_covzm(x, vardim)
-    return cov2cor!(c, collect(sqrt(c[i,i]) for i in 1:min(size(c)...)))
-end
-corzm(x::AbstractVector, y::AbstractMatrix, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt(sum(abs2, x)), sqrt!(sum(abs2, y, dims=vardim)))
-corzm(x::AbstractMatrix, y::AbstractVector, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt(sum(abs2, y)))
-corzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int=1) =
-    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt!(sum(abs2, y, dims=vardim)))
-
-# corm
-
-corm(x::AbstractVector{T}, xmean) where {T} = one(real(T))
-corm(x::AbstractMatrix, xmean, vardim::Int=1) = corzm(x .- xmean, vardim)
-function corm(x::AbstractVector, mx, y::AbstractVector, my)
-    n = length(x)
-    length(y) == n || throw(DimensionMismatch("inconsistent lengths"))
-    n > 0 || throw(ArgumentError("correlation only defined for non-empty vectors"))
-
-    @inbounds begin
-        # Initialize the accumulators
-        xx = zero(sqrt(abs2(x[1])))
-        yy = zero(sqrt(abs2(y[1])))
-        xy = zero(x[1] * y[1]')
-
-        @simd for i in eachindex(x, y)
-            xi = x[i] - mx
-            yi = y[i] - my
-            xx += abs2(xi)
-            yy += abs2(yi)
-            xy += xi * yi'
-        end
-    end
-    return clampcor(xy / max(xx, yy) / sqrt(min(xx, yy) / max(xx, yy)))
-end
-
-corm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean, vardim::Int=1) =
-    corzm(x .- xmean, y .- ymean, vardim)
-
-# cor
-"""
-    cor(x::AbstractVector)
-
-Return the number one.
-"""
-cor(x::AbstractVector) = one(real(eltype(x)))
-
-"""
-    cor(X::AbstractMatrix; dims::Int=1)
-
-Compute the Pearson correlation matrix of the matrix `X` along the dimension `dims`.
-"""
-cor(X::AbstractMatrix; dims::Int=1) = corm(X, _vmean(X, dims), dims)
-
-"""
-    cor(x::AbstractVector, y::AbstractVector)
-
-Compute the Pearson correlation between the vectors `x` and `y`.
-"""
-cor(x::AbstractVector, y::AbstractVector) = corm(x, mean(x), y, mean(y))
-
-"""
-    cor(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims=1)
-
-Compute the Pearson correlation between the vectors or matrices `X` and `Y` along the dimension `dims`.
-"""
-cor(x::AbstractVecOrMat, y::AbstractVecOrMat; dims::Int=1) =
-    corm(x, _vmean(x, dims), y, _vmean(y, dims), dims)
-
-function cov(X::SparseMatrixCSC; dims::Int=1, corrected::Bool=true)
-    vardim = dims
-    a, b = size(X)
-    n, p = vardim == 1 ? (a, b) : (b, a)
-
-    # The covariance can be decomposed into two terms
-    # 1/(n - 1) ∑ (x_i - x̄)*(x_i - x̄)' = 1/(n - 1) (∑ x_i*x_i' - n*x̄*x̄')
-    # which can be evaluated via a sparse matrix-matrix product
-
-    # Compute ∑ x_i*x_i' = X'X using sparse matrix-matrix product
-    out = Matrix(unscaled_covzm(X, vardim))
-
-    # Compute x̄
-    x̄ᵀ = mean(X, dims=vardim)
-
-    # Subtract n*x̄*x̄' from X'X
-    @inbounds for j in 1:p, i in 1:p
-        out[i,j] -= x̄ᵀ[i] * x̄ᵀ[j]' * n
-    end
-
-    # scale with the sample size n or the corrected sample size n - 1
-    return rmul!(out, inv(n - corrected))
-end
-# This is the function that does the reduction underlying var/std
-function centralize_sumabs2!(R::AbstractArray{S}, A::SparseMatrixCSC{Tv,Ti}, means::AbstractArray) where {S,Tv,Ti}
-    lsiz = Base.check_reducedims(R,A)
-    size(means) == size(R) || error("size of means must match size of R")
-    isempty(R) || fill!(R, zero(S))
-    isempty(A) && return R
-
-    colptr = A.colptr
-    rowval = A.rowval
-    nzval = A.nzval
-    m = size(A, 1)
-    n = size(A, 2)
-
-    if size(R, 1) == size(R, 2) == 1
-        # Reduction along both columns and rows
-        R[1, 1] = centralize_sumabs2(A, means[1])
-    elseif size(R, 1) == 1
-        # Reduction along rows
-        @inbounds for col = 1:n
-            mu = means[col]
-            r = convert(S, (m-colptr[col+1]+colptr[col])*abs2(mu))
-            @simd for j = colptr[col]:colptr[col+1]-1
-                r += abs2(nzval[j] - mu)
-            end
-            R[1, col] = r
-        end
-    elseif size(R, 2) == 1
-        # Reduction along columns
-        rownz = fill(convert(Ti, n), m)
-        @inbounds for col = 1:n
-            @simd for j = colptr[col]:colptr[col+1]-1
-                row = rowval[j]
-                R[row, 1] += abs2(nzval[j] - means[row])
-                rownz[row] -= 1
-            end
-        end
-        for i = 1:m
-            R[i, 1] += rownz[i]*abs2(means[i])
-        end
-    else
-        # Reduction along a dimension > 2
-        @inbounds for col = 1:n
-            lastrow = 0
-            @simd for j = colptr[col]:colptr[col+1]-1
-                row = rowval[j]
-                for i = lastrow+1:row-1
-                    R[i, col] = abs2(means[i, col])
-                end
-                R[row, col] = abs2(nzval[j] - means[row, col])
-                lastrow = row
-            end
-            for i = lastrow+1:m
-                R[i, col] = abs2(means[i, col])
-            end
-        end
-    end
-    return R
-end
-
-
-"""
-    linreg(x, y)
-
-Perform simple linear regression using Ordinary Least Squares. Returns `a` and `b` such
-that `a + b*x` is the closest straight line to the given points `(x, y)`, i.e., such that
-the squared error between `y` and `a + b*x` is minimized.
-
-# Examples
-```julia-repl
-julia> x = 1:10
-1:10
-
-julia> y = 2 .* x .+ 5
-7:2:25
-
-julia> a, b = linreg(x, y)
-(5.000000000000002, 1.9999999999999996)
-```
-"""
-function linreg(x::AbstractVector, y::AbstractVector)
-    # Least squares given
-    # Y = a + b*X
-    # where
-    # b = cov(X, Y)/var(X)
-    # a = mean(Y) - b*mean(X)
-    if size(x) != size(y)
-        throw(DimensionMismatch("x has size $(size(x)) and y has size $(size(y)), " *
-            "but these must be the same size"))
-    end
-    mx = mean(x)
-    my = mean(y)
-    # don't need to worry about the scaling (n vs n - 1) since they cancel in the ratio
-    b = covm(x, mx, y, my)/varm(x, mx)
-    a = my - b*mx
-    return (a, b)
-end
diff --git a/src/cov.jl b/src/cov.jl
index 52f30c2f5..773e1bcf2 100644
--- a/src/cov.jl
+++ b/src/cov.jl
@@ -85,7 +85,7 @@ scattermat(x::DenseMatrix, vardim::Int=1) =
     scattermatm(x, mean(x, dims = vardim), vardim)
 
 scattermat(x::DenseMatrix, wv::AbstractWeights, vardim::Int=1) =
-    scattermatm(x, MEANHOME.mean(x, wv, vardim), wv, vardim)
+    scattermatm(x, Statistics.mean(x, wv, vardim), wv, vardim)
 
 ## weighted cov
 covm(x::DenseMatrix, mean, w::AbstractWeights, vardim::Int=1;
diff --git a/src/moments.jl b/src/moments.jl
index e43db4c7f..0d059402b 100644
--- a/src/moments.jl
+++ b/src/moments.jl
@@ -43,7 +43,7 @@ function var(v::RealArray, w::AbstractWeights; mean=nothing,
     corrected = depcheck(:var, corrected)
 
     if mean == nothing
-        varm(v, w, MEANHOME.mean(v, w); corrected=corrected)
+        varm(v, w, Statistics.mean(v, w); corrected=corrected)
     else
         varm(v, w, mean; corrected=corrected)
     end
@@ -66,7 +66,7 @@ function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dim::Int;
         varm!(R, A, w, Base.reducedim_initarray(A, dim, 0, eltype(R)), dim;
                    corrected=corrected)
     elseif mean == nothing
-        varm!(R, A, w, MEANHOME.mean(A, w, dim), dim; corrected=corrected)
+        varm!(R, A, w, Statistics.mean(A, w, dim), dim; corrected=corrected)
     else
         # check size of mean
         for i = 1:ndims(A)
@@ -139,7 +139,7 @@ std(v::RealArray, w::AbstractWeights; mean=nothing, corrected::DepBool=nothing)
     sqrt.(var(v, w; mean=mean, corrected=depcheck(:std, corrected)))
 
 stdm(v::RealArray, m::RealArray, dim::Int; corrected::DepBool=nothing) =
-    sqrt!(StatsCompat.varm(v, m, dims=dim, corrected=depcheck(:stdm, corrected)))
+    sqrt!(varm(v, m, dims=dim, corrected=depcheck(:stdm, corrected)))
 
 stdm(v::RealArray, w::AbstractWeights, m::RealArray, dim::Int;
           corrected::DepBool=nothing) =
@@ -193,7 +193,7 @@ end
 
 function mean_and_var(A::RealArray, dim::Int; corrected::Bool=true)
     m = mean(A, dims = dim)
-    v = StatsCompat.varm(A, m, dims = dim, corrected=corrected)
+    v = varm(A, m, dims = dim, corrected=corrected)
     m, v
 end
 function mean_and_std(A::RealArray, dim::Int; corrected::Bool=true)
diff --git a/src/rankcorr.jl b/src/rankcorr.jl
index 8d326e4d0..ea9077a61 100644
--- a/src/rankcorr.jl
+++ b/src/rankcorr.jl
@@ -20,11 +20,11 @@ of the columns of `x` and `y`.
 corspearman(x::RealVector, y::RealVector) = cor(tiedrank(x), tiedrank(y))
 
 corspearman(X::RealMatrix, Y::RealMatrix) =
-    cor(mapslices(tiedrank, X, 1), mapslices(tiedrank, Y, 1))
-corspearman(X::RealMatrix, y::RealVector) = cor(mapslices(tiedrank, X, 1), tiedrank(y))
-corspearman(x::RealVector, Y::RealMatrix) = cor(tiedrank(x), mapslices(tiedrank, Y, 1))
+    cor(mapslices(tiedrank, X, dims=1), mapslices(tiedrank, Y, dims=1))
+corspearman(X::RealMatrix, y::RealVector) = cor(mapslices(tiedrank, X, dims=1), tiedrank(y))
+corspearman(x::RealVector, Y::RealMatrix) = cor(tiedrank(x), mapslices(tiedrank, Y, dims=1))
 
-corspearman(X::RealMatrix) = (Z = mapslices(tiedrank, X, 1); cor(Z, Z))
+corspearman(X::RealMatrix) = (Z = mapslices(tiedrank, X, dims=1); cor(Z, Z))
 
 
 #######################################
diff --git a/test/ambiguous.jl b/test/ambiguous.jl
index c02f4870c..59014ca0f 100644
--- a/test/ambiguous.jl
+++ b/test/ambiguous.jl
@@ -1,12 +1,8 @@
 using StatsBase
 using Test
-if VERSION >= v"0.7.0-beta.85"
-    using Statistics
-end
+using Statistics
 
 @testset "Ambiguities" begin
-    # Ambiguities with Base and Core introduced by this package
-    tocheck = [StatsBase, Base, Core]
-    VERSION >= v"0.7.0-beta.85" && push!(tocheck, Statistics)
-    @test_broken isempty(detect_ambiguities(tocheck...))
+    # 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/base.jl b/test/base.jl
deleted file mode 100644
index 34a6b8b5b..000000000
--- a/test/base.jl
+++ /dev/null
@@ -1,511 +0,0 @@
-# This file contains code formerly part of Julia. License is MIT: https://julialang.org/license
-
-using Test, Random, LinearAlgebra, SparseArrays
-
-@testset "var & std" begin
-    # edge case: empty vector
-    # iterable; this has to throw for type stability
-    @test_throws ArgumentError var(())
-    @test_throws ArgumentError var((); corrected=false)
-    @test_throws ArgumentError var((); mean=2)
-    @test_throws ArgumentError var((); mean=2, corrected=false)
-    # reduction
-    @test isnan(var(Int[]))
-    @test isnan(var(Int[]; corrected=false))
-    @test isnan(var(Int[]; mean=2))
-    @test isnan(var(Int[]; mean=2, corrected=false))
-    # reduction across dimensions
-    @test isequal(var(Int[], dims=1), [NaN])
-    @test isequal(var(Int[], dims=1; corrected=false), [NaN])
-    @test isequal(var(Int[], dims=1; mean=[2]), [NaN])
-    @test isequal(var(Int[], dims=1; mean=[2], corrected=false), [NaN])
-
-    # edge case: one-element vector
-    # iterable
-    @test isnan(@inferred(var((1,))))
-    @test var((1,); corrected=false) === 0.0
-    @test var((1,); mean=2) === Inf
-    @test var((1,); mean=2, corrected=false) === 1.0
-    # reduction
-    @test isnan(@inferred(var([1])))
-    @test var([1]; corrected=false) === 0.0
-    @test var([1]; mean=2) === Inf
-    @test var([1]; mean=2, corrected=false) === 1.0
-    # reduction across dimensions
-    @test isequal(@inferred(var([1], dims=1)), [NaN])
-    @test var([1], dims=1; corrected=false) ≈ [0.0]
-    @test var([1], dims=1; mean=[2]) ≈ [Inf]
-    @test var([1], dims=1; mean=[2], corrected=false) ≈ [1.0]
-
-    @test var(1:8) == 6.
-    @test varm(1:8,1) == varm(Vector(1:8),1)
-    @test isnan(varm(1:1,1))
-    @test isnan(var(1:1))
-    @test isnan(var(1:-1))
-
-    @test @inferred(var(1.0:8.0)) == 6.
-    @test varm(1.0:8.0,1.0) == varm(Vector(1.0:8.0),1)
-    @test isnan(varm(1.0:1.0,1.0))
-    @test isnan(var(1.0:1.0))
-    @test isnan(var(1.0:-1.0))
-
-    @test @inferred(var(1.0f0:8.0f0)) === 6.f0
-    @test varm(1.0f0:8.0f0,1.0f0) == varm(Vector(1.0f0:8.0f0),1)
-    @test isnan(varm(1.0f0:1.0f0,1.0f0))
-    @test isnan(var(1.0f0:1.0f0))
-    @test isnan(var(1.0f0:-1.0f0))
-
-    @test varm([1,2,3], 2) ≈ 1.
-    @test var([1,2,3]) ≈ 1.
-    @test var([1,2,3]; corrected=false) ≈ 2.0/3
-    @test var([1,2,3]; mean=0) ≈ 7.
-    @test var([1,2,3]; mean=0, corrected=false) ≈ 14.0/3
-
-    @test varm((1,2,3), 2) ≈ 1.
-    @test var((1,2,3)) ≈ 1.
-    @test var((1,2,3); corrected=false) ≈ 2.0/3
-    @test var((1,2,3); mean=0) ≈ 7.
-    @test var((1,2,3); mean=0, corrected=false) ≈ 14.0/3
-    @test_throws ArgumentError var((1,2,3); mean=())
-
-    @test var([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ [2.5 2.5]'
-    @test var([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ [2.0 2.0]'
-
-    @test var(collect(1:99), dims=1) ≈ [825]
-    @test var(Matrix(transpose(collect(1:99))), dims=2) ≈ [825]
-
-    @test stdm([1,2,3], 2) ≈ 1.
-    @test std([1,2,3]) ≈ 1.
-    @test std([1,2,3]; corrected=false) ≈ sqrt(2.0/3)
-    @test std([1,2,3]; mean=0) ≈ sqrt(7.0)
-    @test std([1,2,3]; mean=0, corrected=false) ≈ sqrt(14.0/3)
-
-    @test stdm([1.0,2,3], 2) ≈ 1.
-    @test std([1.0,2,3]) ≈ 1.
-    @test std([1.0,2,3]; corrected=false) ≈ sqrt(2.0/3)
-    @test std([1.0,2,3]; mean=0) ≈ sqrt(7.0)
-    @test std([1.0,2,3]; mean=0, corrected=false) ≈ sqrt(14.0/3)
-
-    @test std([1.0,2,3]; dims=1)[] ≈ 1.
-    @test std([1.0,2,3]; dims=1, corrected=false)[] ≈ sqrt(2.0/3)
-    @test std([1.0,2,3]; dims=1, mean=[0])[] ≈ sqrt(7.0)
-    @test std([1.0,2,3]; dims=1, mean=[0], corrected=false)[] ≈ sqrt(14.0/3)
-
-    @test stdm((1,2,3), 2) ≈ 1.
-    @test std((1,2,3)) ≈ 1.
-    @test std((1,2,3); corrected=false) ≈ sqrt(2.0/3)
-    @test std((1,2,3); mean=0) ≈ sqrt(7.0)
-    @test std((1,2,3); mean=0, corrected=false) ≈ sqrt(14.0/3)
-
-    @test std([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ sqrt.([2.5 2.5]')
-    @test std([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ sqrt.([2.0 2.0]')
-
-    let A = ComplexF64[exp(i*im) for i in 1:10^4]
-        @test varm(A, 0.) ≈ sum(map(abs2, A)) / (length(A) - 1)
-        @test varm(A, mean(A)) ≈ var(A)
-    end
-
-    @test var([1//1, 2//1]) isa Rational{Int}
-    @test var([1//1, 2//1], dims=1) isa Vector{Rational{Int}}
-
-    @test std([1//1, 2//1]) isa Float64
-    @test std([1//1, 2//1], dims=1) isa Vector{Float64}
-end
-
-function safe_cov(x, y, zm::Bool, cr::Bool)
-    n = length(x)
-    if !zm
-        x = x .- mean(x)
-        y = y .- mean(y)
-    end
-    dot(vec(x), vec(y)) / (n - Int(cr))
-end
-X = [1.0  5.0;
-     2.0  4.0;
-     3.0  6.0;
-     4.0  2.0;
-     5.0  1.0]
-Y = [6.0  2.0;
-     1.0  7.0;
-     5.0  8.0;
-     3.0  4.0;
-     2.0  3.0]
-
-@testset "covariance" begin
-    for vd in [1, 2], zm in [true, false], cr in [true, false]
-        # println("vd = $vd: zm = $zm, cr = $cr")
-        if vd == 1
-            k = size(X, 2)
-            Cxx = zeros(k, k)
-            Cxy = zeros(k, k)
-            for i = 1:k, j = 1:k
-                Cxx[i,j] = safe_cov(X[:,i], X[:,j], zm, cr)
-                Cxy[i,j] = safe_cov(X[:,i], Y[:,j], zm, cr)
-            end
-            x1 = vec(X[:,1])
-            y1 = vec(Y[:,1])
-        else
-            k = size(X, 1)
-            Cxx = zeros(k, k)
-            Cxy = zeros(k, k)
-            for i = 1:k, j = 1:k
-                Cxx[i,j] = safe_cov(X[i,:], X[j,:], zm, cr)
-                Cxy[i,j] = safe_cov(X[i,:], Y[j,:], zm, cr)
-            end
-            x1 = vec(X[1,:])
-            y1 = vec(Y[1,:])
-        end
-
-        c = zm ? StatsBase.covm(x1, 0, corrected=cr) :
-                 cov(x1, corrected=cr)
-        @test isa(c, Float64)
-        @test c ≈ Cxx[1,1]
-        @inferred cov(x1, corrected=cr)
-
-        @test cov(X) == StatsBase.covm(X, mean(X, dims=1))
-        C = zm ? StatsBase.covm(X, 0, vd, corrected=cr) :
-                 cov(X, dims=vd, corrected=cr)
-        @test size(C) == (k, k)
-        @test C ≈ Cxx
-        @inferred cov(X, dims=vd, corrected=cr)
-
-        @test cov(x1, y1) == StatsBase.covm(x1, mean(x1), y1, mean(y1))
-        c = zm ? StatsBase.covm(x1, 0, y1, 0, corrected=cr) :
-                 cov(x1, y1, corrected=cr)
-        @test isa(c, Float64)
-        @test c ≈ Cxy[1,1]
-        @inferred cov(x1, y1, corrected=cr)
-
-        if vd == 1
-            @test cov(x1, Y) == StatsBase.covm(x1, mean(x1), Y, mean(Y, dims=1))
-        end
-        C = zm ? StatsBase.covm(x1, 0, Y, 0, vd, corrected=cr) :
-                 cov(x1, Y, dims=vd, corrected=cr)
-        @test size(C) == (1, k)
-        @test vec(C) ≈ Cxy[1,:]
-        @inferred cov(x1, Y, dims=vd, corrected=cr)
-
-        if vd == 1
-            @test cov(X, y1) == StatsBase.covm(X, mean(X, dims=1), y1, mean(y1))
-        end
-        C = zm ? StatsBase.covm(X, 0, y1, 0, vd, corrected=cr) :
-                 cov(X, y1, dims=vd, corrected=cr)
-        @test size(C) == (k, 1)
-        @test vec(C) ≈ Cxy[:,1]
-        @inferred cov(X, y1, dims=vd, corrected=cr)
-
-        @test cov(X, Y) == StatsBase.covm(X, mean(X, dims=1), Y, mean(Y, dims=1))
-        C = zm ? StatsBase.covm(X, 0, Y, 0, vd, corrected=cr) :
-                 cov(X, Y, dims=vd, corrected=cr)
-        @test size(C) == (k, k)
-        @test C ≈ Cxy
-        @inferred cov(X, Y, dims=vd, corrected=cr)
-    end
-end
-
-function safe_cor(x, y, zm::Bool)
-    if !zm
-        x = x .- mean(x)
-        y = y .- mean(y)
-    end
-    x = vec(x)
-    y = vec(y)
-    dot(x, y) / (sqrt(dot(x, x)) * sqrt(dot(y, y)))
-end
-@testset "correlation" begin
-    for vd in [1, 2], zm in [true, false]
-        # println("vd = $vd: zm = $zm")
-        if vd == 1
-            k = size(X, 2)
-            Cxx = zeros(k, k)
-            Cxy = zeros(k, k)
-            for i = 1:k, j = 1:k
-                Cxx[i,j] = safe_cor(X[:,i], X[:,j], zm)
-                Cxy[i,j] = safe_cor(X[:,i], Y[:,j], zm)
-            end
-            x1 = vec(X[:,1])
-            y1 = vec(Y[:,1])
-        else
-            k = size(X, 1)
-            Cxx = zeros(k, k)
-            Cxy = zeros(k, k)
-            for i = 1:k, j = 1:k
-                Cxx[i,j] = safe_cor(X[i,:], X[j,:], zm)
-                Cxy[i,j] = safe_cor(X[i,:], Y[j,:], zm)
-            end
-            x1 = vec(X[1,:])
-            y1 = vec(Y[1,:])
-        end
-
-        c = zm ? StatsBase.corm(x1, 0) : cor(x1)
-        @test isa(c, Float64)
-        @test c ≈ Cxx[1,1]
-        @inferred cor(x1)
-
-        @test cor(X) == StatsBase.corm(X, mean(X, dims=1))
-        C = zm ? StatsBase.corm(X, 0, vd) : cor(X, dims=vd)
-        @test size(C) == (k, k)
-        @test C ≈ Cxx
-        @inferred cor(X, dims=vd)
-
-        @test cor(x1, y1) == StatsBase.corm(x1, mean(x1), y1, mean(y1))
-        c = zm ? StatsBase.corm(x1, 0, y1, 0) : cor(x1, y1)
-        @test isa(c, Float64)
-        @test c ≈ Cxy[1,1]
-        @inferred cor(x1, y1)
-
-        if vd == 1
-            @test cor(x1, Y) == StatsBase.corm(x1, mean(x1), Y, mean(Y, dims=1))
-        end
-        C = zm ? StatsBase.corm(x1, 0, Y, 0, vd) : cor(x1, Y, dims=vd)
-        @test size(C) == (1, k)
-        @test vec(C) ≈ Cxy[1,:]
-        @inferred cor(x1, Y, dims=vd)
-
-        if vd == 1
-            @test cor(X, y1) == StatsBase.corm(X, mean(X, dims=1), y1, mean(y1))
-        end
-        C = zm ? StatsBase.corm(X, 0, y1, 0, vd) : cor(X, y1, dims=vd)
-        @test size(C) == (k, 1)
-        @test vec(C) ≈ Cxy[:,1]
-        @inferred cor(X, y1, dims=vd)
-
-        @test cor(X, Y) == StatsBase.corm(X, mean(X, dims=1), Y, mean(Y, dims=1))
-        C = zm ? StatsBase.corm(X, 0, Y, 0, vd) : cor(X, Y, dims=vd)
-        @test size(C) == (k, k)
-        @test C ≈ Cxy
-        @inferred cor(X, Y, dims=vd)
-    end
-
-    @test cor(repeat(1:17, 1, 17))[2] <= 1.0
-    @test cor(1:17, 1:17) <= 1.0
-    @test cor(1:17, 18:34) <= 1.0
-    let tmp = range(1, stop=85, length=100)
-        tmp2 = Vector(tmp)
-        @test cor(tmp, tmp) <= 1.0
-        @test cor(tmp, tmp2) <= 1.0
-    end
-end
-
-@testset "variance of complex arrays (#13309)" begin
-    z = rand(ComplexF64, 10)
-    @test var(z) ≈ invoke(var, Tuple{Any}, z) ≈ cov(z) ≈ var(z,dims=1)[1] ≈ sum(abs2, z .- mean(z))/9
-    @test isa(var(z), Float64)
-    @test isa(invoke(var, Tuple{Any}, z), Float64)
-    @test isa(cov(z), Float64)
-    @test isa(var(z,dims=1), Vector{Float64})
-    @test varm(z, 0.0) ≈ invoke(varm, Tuple{Any,Float64}, z, 0.0) ≈ sum(abs2, z)/9
-    @test isa(varm(z, 0.0), Float64)
-    @test isa(invoke(varm, Tuple{Any,Float64}, z, 0.0), Float64)
-    @test cor(z) === 1.0
-    v = varm([1.0+2.0im], 0; corrected = false)
-    @test v ≈ 5
-    @test isa(v, Float64)
-end
-
-@testset "cov and cor of complex arrays (issue #21093)" begin
-    x = [2.7 - 3.3im, 0.9 + 5.4im, 0.1 + 0.2im, -1.7 - 5.8im, 1.1 + 1.9im]
-    y = [-1.7 - 1.6im, -0.2 + 6.5im, 0.8 - 10.0im, 9.1 - 3.4im, 2.7 - 5.5im]
-    @test cov(x, y) ≈ 4.8365 - 12.119im
-    @test cov(y, x) ≈ 4.8365 + 12.119im
-    @test cov(x, reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1)
-    @test cov(reshape(x, :, 1), y) ≈ reshape([4.8365 - 12.119im], 1, 1)
-    @test cov(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1)
-    @test cov([x y]) ≈ [21.779 4.8365-12.119im;
-                        4.8365+12.119im 54.548]
-    @test cor(x, y) ≈ 0.14032104449218274 - 0.35160772008699703im
-    @test cor(y, x) ≈ 0.14032104449218274 + 0.35160772008699703im
-    @test cor(x, reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
-    @test cor(reshape(x, :, 1), y) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
-    @test cor(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1)
-    @test cor([x y]) ≈ [1.0                                          0.14032104449218274-0.35160772008699703im
-                        0.14032104449218274+0.35160772008699703im  1.0]
-end
-
-@testset "Issue #17153 and PR #17154" begin
-    a = rand(10,10)
-    b = copy(a)
-    x = var(a, dims=1)
-    @test b == a
-    x = var(a, dims=2)
-    @test b == a
-    x = std(a, dims=1)
-    @test b == a
-    x = std(a, dims=2)
-    @test b == a
-end
-
-@testset "Unitful elements" begin
-    r = Furlong(1):Furlong(1):Furlong(2)
-    a = Vector(r)
-    @test var(r) == var(a) == Furlong{2}(0.5)
-    @test std(r) == std(a) == Furlong{1}(sqrt(0.5))
-
-    # Issue #21786
-    A = [Furlong{1}(rand(-5:5)) for i in 1:2, j in 1:2]
-    @test var(A, dims=1)[1] === var(A[:, 1])
-    @test std(A, dims=1)[1] === std(A[:, 1])
-end
-
-# Issue #22901
-@testset "var and std of Any arrays" begin
-    x = Any[1, 2, 4, 10]
-    y = Any[1, 2, 4, 10//1]
-    @test var(x) === 16.25
-    @test var(y) === 65//4
-    @test std(x) === sqrt(16.25)
-end
-
-@testset "Promotion in covzm. Issue #8080" begin
-    A = [1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
-    @test StatsBase.covzm(A) - mean(A, dims=1)'*mean(A, dims=1)*size(A, 1)/(size(A, 1) - 1) ≈ cov(A)
-    A = [1//1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1]
-    @test (A'A - size(A, 1)*Base.mean(A, dims=1)'*Base.mean(A, dims=1))/4 == cov(A)
-end
-
-@testset "cov/var/std of Vector{Vector}" begin
-    x = [[2,4,6],[4,6,8]]
-    @test var(x) ≈ vec(var([x[1] x[2]], dims=2))
-    @test std(x) ≈ vec(std([x[1] x[2]], dims=2))
-    @test cov(x) ≈ cov([x[1] x[2]], dims=2)
-end
-
-# Faster covariance function for sparse matrices
-# Prevents densifying the input matrix when subtracting the mean
-# Test against dense implementation
-# PR https://github.com/JuliaLang/julia/pull/22735
-# Part of this test needed to be hacked due to the treatment
-# of Inf in sparse matrix algebra
-# https://github.com/JuliaLang/julia/issues/22921
-# The issue will be resolved in
-# https://github.com/JuliaLang/julia/issues/22733
-@testset "optimizing sparse $elty covariance" for elty in (Float64, Complex{Float64})
-    n = 10
-    p = 5
-    np2 = div(n*p, 2)
-    nzvals, x_sparse = guardsrand(1) do
-        if elty <: Real
-            nzvals = randn(np2)
-        else
-            nzvals = complex.(randn(np2), randn(np2))
-        end
-        nzvals, sparse(rand(1:n, np2), rand(1:p, np2), nzvals, n, p)
-    end
-    x_dense  = convert(Matrix{elty}, x_sparse)
-    @testset "Test with no Infs and NaNs, vardim=$vardim, corrected=$corrected" for vardim in (1, 2),
-                                                                                 corrected in (true, false)
-        @test cov(x_sparse, dims=vardim, corrected=corrected) ≈
-              cov(x_dense , dims=vardim, corrected=corrected)
-    end
-
-    @testset "Test with $x11, vardim=$vardim, corrected=$corrected" for x11 in (NaN, Inf),
-                                                                     vardim in (1, 2),
-                                                                  corrected in (true, false)
-        x_sparse[1,1] = x11
-        x_dense[1 ,1] = x11
-
-        cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected)
-        cov_dense  = cov(x_dense , dims=vardim, corrected=corrected)
-        @test cov_sparse[2:end, 2:end] ≈ cov_dense[2:end, 2:end]
-        @test isfinite.(cov_sparse) == isfinite.(cov_dense)
-        @test isfinite.(cov_sparse) == isfinite.(cov_dense)
-    end
-
-    @testset "Test with NaN and Inf, vardim=$vardim, corrected=$corrected" for vardim in (1, 2),
-                                                                            corrected in (true, false)
-        x_sparse[1,1] = Inf
-        x_dense[1 ,1] = Inf
-        x_sparse[2,1] = NaN
-        x_dense[2 ,1] = NaN
-
-        cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected)
-        cov_dense  = cov(x_dense , dims=vardim, corrected=corrected)
-        @test cov_sparse[(1 + vardim):end, (1 + vardim):end] ≈
-              cov_dense[ (1 + vardim):end, (1 + vardim):end]
-        @test isfinite.(cov_sparse) == isfinite.(cov_dense)
-        @test isfinite.(cov_sparse) == isfinite.(cov_dense)
-    end
-end
-
-# issue #6672
-@test std(AbstractFloat[1,2,3], dims=1) == [1.0]
-
-@testset "var of sparse array" begin
-    se33 = SparseMatrixCSC{Float64}(I, 3, 3)
-    sA = sprandn(3, 7, 0.5)
-    pA = sparse(rand(3, 7))
-
-    for arr in (se33, sA, pA)
-        farr = Array(arr)
-        @test var(arr) ≈ var(farr)
-        @test var(arr, dims=1) ≈ var(farr, dims=1)
-        @test var(arr, dims=2) ≈ var(farr, dims=2)
-        @test var(arr, dims=(1, 2)) ≈ [var(farr)]
-        @test isequal(var(arr, dims=3), var(farr, dims=3))
-    end
-
-    @testset "empty cases" begin
-        @test var(sparse(Int[])) === NaN
-        @test isequal(var(spzeros(0, 1), dims=1), var(Matrix{Int}(I, 0, 1), dims=1))
-        @test isequal(var(spzeros(0, 1), dims=2), var(Matrix{Int}(I, 0, 1), dims=2))
-        @test isequal(var(spzeros(0, 1), dims=(1, 2)), var(Matrix{Int}(I, 0, 1), dims=(1, 2)))
-        @test isequal(var(spzeros(0, 1), dims=3), var(Matrix{Int}(I, 0, 1), dims=3))
-    end
-end
-
-@testset "var: empty cases" begin
-    A = Matrix{Int}(undef, 0,1)
-    @test var(A) === NaN
-
-    @test isequal(var(A, dims=1), fill(NaN, 1, 1))
-    @test isequal(var(A, dims=2), fill(NaN, 0, 1))
-    @test isequal(var(A, dims=(1, 2)), fill(NaN, 1, 1))
-    @test isequal(var(A, dims=3), fill(NaN, 0, 1))
-end
-
-@testset "linreg" begin
-    # make sure unequal input arrays throw an error
-    x = [2; 5; 6]
-    y = [3; 7; 10; 10]
-    @test_throws DimensionMismatch linreg(x, y)
-    x = [2 5 6]
-    y = [3; 7; 10]
-    @test_throws MethodError linreg(x, y)
-
-    # check (UnitRange, Array)
-    x = 1:12
-    y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99]
-    @test [linreg(x,y)...] ≈ [2.5559090909090867, 1.6960139860139862]
-    @test [linreg(view(x,1:6),view(y,1:6))...] ≈ [3.8366666666666642,1.3271428571428574]
-
-    # check (LinRange, UnitRange)
-    x = range(1.0, stop=12.0, length=100)
-    y = -100:-1
-    @test [linreg(x, y)...] ≈ [-109.0, 9.0]
-
-    # check (UnitRange, UnitRange)
-    x = 1:12
-    y = 12:-1:1
-    @test [linreg(x, y)...] ≈ [13.0, -1.0]
-
-    # check (LinRange, LinRange)
-    x = range(-5, stop=10, length=100)
-    y = range(50, stop=200, length=100)
-    @test [linreg(x, y)...] ≈ [100.0, 10.0]
-
-    # check (Array, Array)
-    # Anscombe's quartet (https://en.wikipedia.org/wiki/Anscombe%27s_quartet)
-    x123 = [10.0; 8.0; 13.0; 9.0; 11.0; 14.0; 6.0; 4.0; 12.0; 7.0; 5.0]
-    y1 = [8.04; 6.95; 7.58; 8.81; 8.33; 9.96; 7.24; 4.26; 10.84; 4.82; 5.68]
-    @test [linreg(x123,y1)...] ≈ [3.0,0.5] atol=15e-5
-
-    y2 = [9.14; 8.14; 8.74; 8.77; 9.26; 8.10; 6.12; 3.10; 9.13; 7.26; 4.74]
-    @test [linreg(x123,y2)...] ≈ [3.0,0.5] atol=10e-3
-
-    y3 = [7.46; 6.77; 12.74; 7.11; 7.81; 8.84; 6.08; 5.39; 8.15; 6.42; 5.73]
-    @test [linreg(x123,y3)...] ≈ [3.0,0.5] atol=10e-3
-
-    x4 = [8.0; 8.0; 8.0; 8.0; 8.0; 8.0; 8.0; 19.0; 8.0; 8.0; 8.0]
-    y4 = [6.58; 5.76; 7.71; 8.84; 8.47; 7.04; 5.25; 12.50; 5.56; 7.91; 6.89]
-    @test [linreg(x4,y4)...] ≈ [3.0,0.5] atol=10e-3
-end
diff --git a/test/cov.jl b/test/cov.jl
index f96268cc9..fbe560e40 100644
--- a/test/cov.jl
+++ b/test/cov.jl
@@ -85,15 +85,15 @@ weight_funcs = (weights, aweights, fweights, pweights)
         @testset "Mean and covariance" begin
             (m, C) = mean_and_cov(X; corrected=false)
             @test m == mean(X, dims = 1)
-            @test C == StatsCompat.cov(X, dims = 1, corrected=false)
+            @test C == cov(X, dims = 1, corrected=false)
 
             (m, C) = mean_and_cov(X, 1; corrected=false)
             @test m == mean(X, dims = 1)
-            @test C == StatsCompat.cov(X, dims = 1, corrected = false)
+            @test C == cov(X, dims = 1, corrected = false)
 
             (m, C) = mean_and_cov(X, 2; corrected=false)
             @test m == mean(X, dims = 2)
-            @test C == StatsCompat.cov(X, dims = 2, corrected = false)
+            @test C == cov(X, dims = 2, corrected = false)
 
             (m, C) = mean_and_cov(X, wv1; corrected=false)
             @test m == mean(X, wv1, 1)
@@ -118,14 +118,14 @@ weight_funcs = (weights, aweights, fweights, pweights)
             cor2 = cor(X, wv2, 2)
 
             @testset "cov2cor" begin
-                @test cov2cor(StatsCompat.cov(X, dims = 1), StatsCompat.std(X, dims = 1)) ≈ StatsCompat.cor(X, dims = 1)
-                @test cov2cor(StatsCompat.cov(X, dims = 2), StatsCompat.std(X, dims = 2)) ≈ StatsCompat.cor(X, dims = 2)
+                @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
+                @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
                 @test cov2cor(cov1, std1) ≈ cor1
                 @test cov2cor(cov2, std2) ≈ cor2
             end
             @testset "cor2cov" begin
-                @test cor2cov(StatsCompat.cor(X, dims = 1), StatsCompat.std(X, dims = 1)) ≈ StatsCompat.cov(X, dims = 1)
-                @test cor2cov(StatsCompat.cor(X, dims = 2), StatsCompat.std(X, dims = 2)) ≈ StatsCompat.cov(X, dims = 2)
+                @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
+                @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
                 @test cor2cov(cor1, std1) ≈ cov1
                 @test cor2cov(cor2, std2) ≈ cov2
             end
@@ -156,15 +156,15 @@ weight_funcs = (weights, aweights, fweights, pweights)
         @testset "Mean and covariance" begin
             (m, C) = mean_and_cov(X; corrected=true)
             @test m == mean(X, dims =1)
-            @test C == StatsCompat.cov(X, dims = 1, corrected = true)
+            @test C == cov(X, dims = 1, corrected = true)
 
             (m, C) = mean_and_cov(X, 1; corrected=true)
             @test m == mean(X, dims = 1)
-            @test C == StatsCompat.cov(X, dims = 1, corrected = true)
+            @test C == cov(X, dims = 1, corrected = true)
 
             (m, C) = mean_and_cov(X, 2; corrected=true)
             @test m == mean(X, dims = 2)
-            @test C == StatsCompat.cov(X, dims = 2, corrected = true)
+            @test C == cov(X, dims = 2, corrected = true)
 
             if isa(wv1, Weights)
                 @test_throws ArgumentError mean_and_cov(X, wv1; corrected=true)
@@ -194,8 +194,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
                 cor2 = cor(X, wv2, 2)
 
                 @testset "cov2cor" begin
-                    @test cov2cor(StatsCompat.cov(X, dims = 1), StatsCompat.std(X, dims = 1)) ≈ StatsCompat.cor(X, dims = 1)
-                    @test cov2cor(StatsCompat.cov(X, dims = 2), StatsCompat.std(X, dims = 2)) ≈ StatsCompat.cor(X, dims = 2)
+                    @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
+                    @test cov2cor(cov(X, dims = 2), std(X, dims = 2)) ≈ cor(X, dims = 2)
                     @test cov2cor(cov1, std1) ≈ cor1
                     @test cov2cor(cov2, std2) ≈ cor2
                 end
@@ -213,8 +213,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
                 end
 
                 @testset "cor2cov" begin
-                    @test cor2cov(StatsCompat.cor(X, dims = 1), StatsCompat.std(X, dims = 1)) ≈ StatsCompat.cov(X, dims = 1)
-                    @test cor2cov(StatsCompat.cor(X, dims = 2), StatsCompat.std(X, dims = 2)) ≈ StatsCompat.cov(X, dims = 2)
+                    @test cor2cov(cor(X, dims = 1), std(X, dims = 1)) ≈ cov(X, dims = 1)
+                    @test cor2cov(cor(X, dims = 2), std(X, dims = 2)) ≈ cov(X, dims = 2)
                     @test cor2cov(cor1, std1) ≈ cov1
                     @test cor2cov(cor2, std2) ≈ cov2
                 end
@@ -235,8 +235,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
     end
 
     @testset "Correlation" begin
-        @test cor(X, f(ones(3)), 1) ≈ StatsCompat.cor(X, dims = 1)
-        @test cor(X, f(ones(8)), 2) ≈ StatsCompat.cor(X, dims = 2)
+        @test cor(X, f(ones(3)), 1) ≈ cor(X, dims = 1)
+        @test cor(X, f(ones(8)), 2) ≈ cor(X, dims = 2)
 
         cov1 = cov(X, wv1, 1; corrected=false)
         std1 = std(X, wv1, 1; corrected=false)
diff --git a/test/dimensionful.jl b/test/dimensionful.jl
deleted file mode 100644
index 5ac003727..000000000
--- a/test/dimensionful.jl
+++ /dev/null
@@ -1,74 +0,0 @@
-# This file was formerly part of Julia. License is MIT: https://julialang.org/license
-
-# Here we implement a minimal dimensionful type Furlong, which is used
-# to test dimensional correctness of various functions in Base.  Furlong
-# is exported by the TestHelpers module.
-
-# represents a quantity in furlongs^p
-struct Furlong{p,T<:Number} <: Number
-    val::T
-    Furlong{p,T}(v::Number) where {p,T} = new(v)
-    Furlong{p,T}(x::Furlong{p}) where {p,T} = new(x.val)
-end
-Furlong(x::T) where {T<:Number} = Furlong{1,T}(x)
-(::Type{T})(x::Furlong) where {T<:Number} = T(x.val)
-Furlong{p}(v::Number) where {p} = Furlong{p,typeof(v)}(v)
-Furlong{p,T}(x::Furlong{p,S}) where {T,p,S} = Furlong{p,T}(T(x.val))
-
-Base.promote_type(::Type{Furlong{p,T}}, ::Type{Furlong{p,S}}) where {p,T,S} =
-    (Base.@_pure_meta; Furlong{p,promote_type(T,S)})
-
-Base.one(x::Furlong{p,T}) where {p,T} = one(T)
-Base.one(::Type{Furlong{p,T}}) where {p,T} = one(T)
-Base.zero(x::Furlong{p,T}) where {p,T} = Furlong{p,T}(zero(T))
-Base.zero(::Type{Furlong{p,T}}) where {p,T} = Furlong{p,T}(zero(T))
-Base.iszero(x::Furlong) = iszero(x.val)
-
-# convert Furlong exponent p to a canonical form.  This
-# is not type stable, but it doesn't matter since it is used
-# at compile time (in generated functions), not runtime
-canonical_p(p) = isinteger(p) ? Int(p) : Rational{Int}(p)
-
-Base.abs(x::Furlong{p}) where {p} = Furlong{p}(abs(x.val))
-@generated Base.abs2(x::Furlong{p}) where {p} = :(Furlong{$(canonical_p(2p))}(abs2(x.val)))
-@generated Base.inv(x::Furlong{p}) where {p} = :(Furlong{$(canonical_p(-p))}(inv(x.val)))
-LinearAlgebra.sylvester(a::Furlong,b::Furlong,c::Furlong) = -c / (a + b)
-
-for f in (:isfinite, :isnan, :isreal)
-    @eval Base.$f(x::Furlong) = $f(x.val)
-end
-for f in (:real,:imag,:complex,:middle,:+,:-)
-    @eval Base.$f(x::Furlong{p}) where {p} = Furlong{p}($f(x.val))
-end
-
-for op in (:+, :-, :middle, :hypot)
-    @eval function Base.$op(x::Furlong{p}, y::Furlong{p}) where {p}
-        v = $op(x.val, y.val)
-        Furlong{p}(v)
-    end
-end
-for op in (:(==), :(!=), :<, :<=, :isless, :isequal)
-    @eval Base.$op(x::Furlong{p}, y::Furlong{p}) where {p} = $op(x.val, y.val)
-end
-# generated functions to allow type inference of the value of the exponent:
-for (f,op) in ((:_plus,:+),(:_minus,:-),(:_times,:*),(:_div,://))
-    @eval @generated function $f(v::T, ::Furlong{p}, ::Union{Furlong{q},Val{q}}) where {T,p,q}
-        s = $op(p, q)
-        :(Furlong{$(canonical_p(s)),$T}(v))
-    end
-end
-for (op,eop) in ((:*, :_plus), (:/, :_minus), (://, :_minus), (:div, :_minus))
-    @eval begin
-        Base.$op(x::Furlong{p}, y::Furlong{q}) where {p,q} =
-            $eop($op(x.val, y.val),x,y)
-        Base.$op(x::Furlong{p}, y::S) where {p,S<:Number} = $op(x,Furlong{0,S}(y))
-        Base.$op(x::S, y::Furlong{p}) where {p,S<:Number} = $op(Furlong{0,S}(x),y)
-    end
-end
-for op in (:rem, :mod)
-    @eval begin
-        Base.$op(x::Furlong{p}, y::Furlong) where {p} = Furlong{p}($op(x.val, y.val))
-        Base.$op(x::Furlong{p}, y::Number) where {p} = Furlong{p}($op(x.val, y))
-    end
-end
-Base.sqrt(x::Furlong) = _div(sqrt(x.val), x, Val(2))
diff --git a/test/moments.jl b/test/moments.jl
index 9e03b9185..dbf812059 100644
--- a/test/moments.jl
+++ b/test/moments.jl
@@ -136,7 +136,7 @@ w2 = rand(6)
         for d in 1:2
             (m, v) = mean_and_var(x, d; corrected=false)
             @test m == mean(x, dims = d)
-            @test v == StatsCompat.var(x, dims = d, corrected=false)
+            @test v == var(x, dims = d, corrected=false)
         end
 
         (m, v) = mean_and_var(x, wv1, 1; corrected=false)
@@ -152,7 +152,7 @@ w2 = rand(6)
         for d in 1:2
             (m, s) = mean_and_std(x, d; corrected=false)
             @test m == mean(x, dims = d)
-            @test s == StatsCompat.std(x, dims = d; corrected=false)
+            @test s == std(x, dims = d; corrected=false)
         end
 
         (m, s) = mean_and_std(x, wv1, 1; corrected=false)
@@ -204,7 +204,7 @@ end
         for d in 1:2
             (m, v) = mean_and_var(x, d; corrected=true)
             @test m == mean(x, dims = d)
-            @test v == StatsCompat.var(x, dims = d, corrected=true)
+            @test v == var(x, dims = d, corrected=true)
         end
 
         if isa(wv1, Weights)
@@ -224,7 +224,7 @@ end
         for d in 1:2
             (m, s) = mean_and_std(x, d; corrected=true)
             @test m == mean(x, dims = d)
-            @test s == StatsCompat.std(x, dims = d, corrected=true)
+            @test s == std(x, dims = d, corrected=true)
         end
 
         if isa(wv1, Weights)
diff --git a/test/runtests.jl b/test/runtests.jl
index e05781067..cbc5e6dfa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,9 +1,7 @@
 using StatsBase
 using LinearAlgebra
 using Random
-if VERSION >= v"0.7.0-beta.85"
-    using Statistics
-end
+using Statistics
 
 tests = ["ambiguous",
          "weights",
@@ -21,15 +19,9 @@ tests = ["ambiguous",
          "robust",
          "sampling",
          "wsampling",
-         "statmodels",
-         "statscompat"]#,
+         "statmodels"]
          #"statquiz"]
 
-if StatsBase.BASESTATS_IN_STATSBASE
-    push!(tests, "base")
-    include("dimensionful.jl")
-end
-
 println("Running tests:")
 
 for t in tests
diff --git a/test/sampling.jl b/test/sampling.jl
index d37cbcf85..10c5c374e 100644
--- a/test/sampling.jl
+++ b/test/sampling.jl
@@ -78,19 +78,11 @@ test_rng_use(sample, 1:10, 10)
 
     srand(1)
 
-    if VERSION < v"0.7.0-DEV.1325" # new faster rand(MersenneTwister, UnitRange)
-        @test samplepair(2)  ===  (1, 2)
-        @test samplepair(10) === (10, 6)
+    @test samplepair(2)  ===  (1, 2)
+    @test samplepair(10) === (8, 2)
 
-        @test samplepair([3, 4, 2, 6, 8]) === (6, 2)
-        @test samplepair([1, 2])          === (1, 2)
-    else
-        @test samplepair(2)  ===  (1, 2)
-        @test samplepair(10) === (8, 2)
-
-        @test samplepair([3, 4, 2, 6, 8]) === (2, 6)
-        @test samplepair([1, 2])          === (1, 2)
-    end
+    @test samplepair([3, 4, 2, 6, 8]) === (2, 6)
+    @test samplepair([1, 2])          === (1, 2)
 end
 
 test_rng_use(samplepair, 1000)
diff --git a/test/scalarstats.jl b/test/scalarstats.jl
index 21973bda6..9d7bf3a5e 100644
--- a/test/scalarstats.jl
+++ b/test/scalarstats.jl
@@ -58,8 +58,8 @@ z2 = [8. 2. 3. 1.; 24. 10. -1. -1.; 20. 12. 1. -2.]
 @test zscore!(zeros(size(a)), a, [1 3 2 4], [0.25 0.5 1.0 2.0]) ≈ z2
 
 @test zscore(a)    ≈ zscore(a, mean(a), std(a))
-@test zscore(a, 1) ≈ zscore(a, mean(a, dims=1), StatsCompat.std(a, dims=1))
-@test zscore(a, 2) ≈ zscore(a, mean(a, dims=2), StatsCompat.std(a, dims=2))
+@test zscore(a, 1) ≈ zscore(a, mean(a, dims=1), std(a, dims=1))
+@test zscore(a, 2) ≈ zscore(a, mean(a, dims=2), std(a, dims=2))
 
 
 ###### quantile & friends
diff --git a/test/statscompat.jl b/test/statscompat.jl
deleted file mode 100644
index d71d73ec2..000000000
--- a/test/statscompat.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-using StatsBase, Test
-
-@testset "StatsCompat for replacing some Compat functions" begin
-    @test StatsCompat.varm([1 2; 3 4], -1) == 18
-    @test StatsCompat.varm([1 2; 3 4], [-1 -2], dims=1) == [20 52]
-    @test StatsCompat.varm([1 2; 3 4], [-1, -2], dims=2) == hcat([13, 61])
-    @test StatsCompat.var([1 2; 3 4]) == 5/3
-    @test StatsCompat.var([1 2; 3 4], dims=1) == [2 2]
-    @test StatsCompat.var([1 2; 3 4], dims=2) == hcat([0.5, 0.5])
-    @test StatsCompat.var([1 2; 3 4], corrected=false) == 1.25
-    @test StatsCompat.var([1 2; 3 4], corrected=false, dims=1) == [1 1]
-    @test StatsCompat.var([1 2; 3 4], corrected=false, dims=2) == hcat([0.25, 0.25])
-    @test StatsCompat.std([1 2; 3 4]) == sqrt(5/3)
-    @test StatsCompat.std([1 2; 3 4], dims=1) == [sqrt(2) sqrt(2)]
-    @test StatsCompat.std([1 2; 3 4], dims=2) == hcat([sqrt(0.5), sqrt(0.5)])
-    @test StatsCompat.std([1 2; 3 4], corrected=false) == sqrt(1.25)
-    @test StatsCompat.std([1 2; 3 4], corrected=false, dims=1) == [sqrt(1) sqrt(1)]
-    @test StatsCompat.std([1 2; 3 4], corrected=false, dims=2) == hcat([sqrt(0.25), sqrt(0.25)])
-    @test StatsCompat.cov([1 2; 3 4]) == [2 2; 2 2]
-    @test StatsCompat.cov([1 2; 3 4], dims=1) == [2 2; 2 2]
-    @test StatsCompat.cov([1 2; 3 4], dims=2) == [0.5 0.5; 0.5 0.5]
-    @test StatsCompat.cov([1 2; 3 4], [4; 5]) == hcat([1; 1])
-    @test StatsCompat.cov([1 2; 3 4], [4; 5], dims=1) == hcat([1; 1])
-    @test StatsCompat.cov([1 2; 3 4], [4; 5], dims=2) == hcat([0.5; 0.5])
-    @test StatsCompat.cov([1 2; 3 4], [4; 5], corrected=false) == hcat([0.5; 0.5])
-    @test StatsCompat.cov([1 2; 3 4], [4; 5], corrected=false, dims=1) == hcat([0.5; 0.5])
-    @test StatsCompat.cov([1 2; 3 4], [4; 5], corrected=false, dims=2) == hcat([0.25; 0.25])
-    @test StatsCompat.cor([1 2; 3 4]) ≈ [1 1; 1 1]
-    @test StatsCompat.cor([1 2; 3 4], dims=1) ≈ [1 1; 1 1]
-    @test StatsCompat.cor([1 2; 3 4], dims=2) ≈ [1 1; 1 1]
-    @test StatsCompat.cor([1 2; 3 4], [4; 5]) ≈ [1; 1]
-    @test StatsCompat.cor([1 2; 3 4], [4; 5], dims=1) ≈ [1; 1]
-    @test StatsCompat.cor([1 2; 3 4], [4; 5], dims=2) ≈ [1; 1]
-end