diff --git a/Project.toml b/Project.toml
index d5043c2..55a57f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,8 @@ authors = ["Phillip Alday <me@phillipalday.com> and contributors"]
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
@@ -20,6 +22,8 @@ GLM = "1.9"
 LinearAlgebra = "1"
 MixedModels = "5.9"
 MixedModelsDatasets = "0.2"
+ProgressMeter = "1"
+Random = "1"
 RDatasets = "0.7.7, 0.8"
 StableRNGs = "1"
 Statistics = "1"
diff --git a/src/MixedModelsExtras.jl b/src/MixedModelsExtras.jl
index a045a5c..c1bffbf 100644
--- a/src/MixedModelsExtras.jl
+++ b/src/MixedModelsExtras.jl
@@ -2,11 +2,13 @@ module MixedModelsExtras
 
 using LinearAlgebra
 using MixedModels
+using Random
 using Statistics
 using StatsBase
 using StatsModels
 using Tables
 
+using MixedModels: replicate
 using GLM: linkinv, Link
 using StatsModels: termnames, vif, gvif
 export termnames, gvif, vif
@@ -25,6 +27,9 @@ export partial_fitted
 include("shrinkage.jl")
 export shrinkagenorm, shrinkagetables
 
+include("bootstrap.jl")
+export bootstrap_lrt
+
 include("tables.jl")
 export ictable
 
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
new file mode 100644
index 0000000..7f2d0ff
--- /dev/null
+++ b/src/bootstrap.jl
@@ -0,0 +1,84 @@
+"""
+    bootstrap_lrt(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
+                  optsum_overrides=(;), progress=true)
+
+Bootstrapped likelihood ratio test applied to a set of nested models.
+
+The first model is used to simulate `n` dataset replicates, where the ground truth is that
+specified by the first model. Each of the other models is then refit to those
+null data and the underlying distribution of deviance differences is then captured.
+For final computation of the p-values, the observed difference in deviance  between the
+original models is compared against this null distribution.
+
+!!! note
+    The precision of the resulting p-value cannot exceed ``1/n``.
+
+!!! warn
+    This method is **not** thread safe. For efficiency , the models are modified
+    during bootstrapping and the original fits are only restored at the end.
+
+!!! note
+    The nesting of the models is not checked.  It is incumbent on the user
+    to check this. This differs from `StatsModels.lrtest` as nesting in
+    mixed models, especially in the random effects specification, may be non obvious.
+
+This functionality may be deprecated in the future in favor of `StatsModels.lrtest`.
+"""
+function bootstrap_lrt(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
+                       optsum_overrides=(;), progress=true)
+    y0 = copy(response(m0))
+    ys = [copy(response(m)) for m in ms]
+    local models
+    local dofs
+    local formulas
+    local devs
+    local lls
+    local devdiffs
+    local pvals
+    try
+        models = [m0; ms...]
+        models = models[sortperm(dof.(models))]
+        dofs = dof.(models)
+        formulas = string.(formula.(models))
+        devs = deviance.(models)
+        lls = loglikelihood.(models)
+        devdiffs = .-(diff(devs))
+
+        for (key, val) in pairs(optsum_overrides)
+            setfield!(m0.optsum, key, val)
+            for m in ms
+                setfield!(m.optsum, key, val)
+            end
+        end
+        nulldist = replicate(n; progress) do
+            simulate!(rng, m0)
+            refit!(m0; progress=false)
+            for m in ms
+                refit!(m, response(m0); progress=false)
+            end
+            return deviance.(models)
+        end
+        nulldist = stack(nulldist; dims=1)
+        nulldist = .-(diff(nulldist; dims=2))
+        pvals = map(enumerate(devdiffs)) do (idx, dev)
+            if dev > 0
+                mean(>=(dev), view(nulldist, :, idx))
+            else
+                NaN
+            end
+        end
+    finally
+        # restore the original fits
+        if progress
+            @info "Bootstrapping complete, cleaning up..."
+        end
+        refit!(m0, y0; progress=false)
+        for (m, y) in zip(ms, ys)
+            refit!(m, y; progress=false)
+        end
+    end
+    lrt = StatsModels.LRTestResult(Int(nobs(m0)), Tuple(devs), Tuple(lls), Tuple(dofs),
+                                   (NaN, pvals...))
+    return MixedModels.LikelihoodRatioTest(Tuple(formulas), lrt,
+                                           first(models) isa LinearMixedModel)
+end
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
new file mode 100644
index 0000000..d53221f
--- /dev/null
+++ b/test/bootstrap.jl
@@ -0,0 +1,46 @@
+sleepstudy = dataset(:sleepstudy)
+fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | subj)),
+          sleepstudy; progress)
+fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)),
+          sleepstudy; progress)
+fmzc = fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days | subj)),
+           sleepstudy; progress)
+
+# fm0 ⊂ fmzc ⊂ fm1, in ascending dof order, so the analytic LRT is well defined
+# and gives us a reference to check the bootstrap version against.
+lrt = likelihoodratiotest(fm0, fmzc, fm1)
+d0, dzc, d1 = deviance(fm0), deviance(fmzc), deviance(fm1)
+
+boot = bootstrap_lrt(StableRNG(42), 200, fm0, fmzc, fm1;
+                     progress, optsum_overrides=(; maxfeval=500))
+
+@testset "bookkeeping matches the analytic LRT" begin
+    @test boot.dof == lrt.dof
+    @test all(isapprox.(boot.deviance, lrt.deviance))
+    @test boot.formulas == lrt.formulas
+    @test boot.linear
+end
+
+@testset "p-values track the sign of the analytic LRT" begin
+    @test isnan(boot.pvalues[1])
+    # fm0 vs fmzc: huge deviance drop, analytic p ≈ 9e-11
+    @test boot.pvalues[2] < 0.05
+    # fmzc vs fm1: negligible deviance drop, analytic p ≈ 0.80
+    @test boot.pvalues[3] > 0.3
+end
+
+@testset "original fits are restored after bootstrapping" begin
+    @test deviance(fm0) ≈ d0
+    @test deviance(fmzc) ≈ dzc
+    @test deviance(fm1) ≈ d1
+end
+
+@testset "bookkeeping is independent of the order of ms..." begin
+    # only the identity of m0 (the data-generating model) should matter;
+    # re-sorting by dof for reporting must not depend on argument order
+    boot2 = bootstrap_lrt(StableRNG(42), 200, fm0, fm1, fmzc; progress)
+    @test boot2.dof == boot.dof
+    @test boot2.formulas == boot.formulas
+    @test boot2.deviance == boot.deviance
+    @test isequal(boot2.pvalues, boot.pvalues)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index d3a1c90..8b31bc2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -28,6 +28,10 @@ end
     include("shrinkage.jl")
 end
 
+@testset "bootstrap" begin
+    include("bootstrap.jl")
+end
+
 @testset "tables" begin
     include("tables.jl")
 end
diff --git a/test/set_up_tests.jl b/test/set_up_tests.jl
index 494c3e9..79af18f 100644
--- a/test/set_up_tests.jl
+++ b/test/set_up_tests.jl
@@ -13,6 +13,7 @@ using Tables
 using Test
 
 using GLM: linkinv, Link
+using MixedModels: likelihoodratiotest
 using MixedModelsDatasets: dataset
 using MixedModelsExtras: _ranef
 using RDatasets: dataset as rdataset