From 57544873677036ff24e31520b20e3ac8bef02448 Mon Sep 17 00:00:00 2001
From: Phillip Alday <phillip.alday@beacon.bio>
Date: Tue, 19 Apr 2022 10:48:57 -0500
Subject: [PATCH 01/14] icc tolerances

---
 test/icc.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/icc.jl b/test/icc.jl
index 163fcc5..9bccacf 100644
--- a/test/icc.jl
+++ b/test/icc.jl
@@ -11,7 +11,8 @@ progress = false
     model = fit(MixedModel, @formula(reaction ~ 1 + (1 | subj)), dataset(:sleepstudy);
                 progress)
     @test icc(model, :subj) == icc(model, [:subj]) == icc(model)
-    @test icc(model, :subj) ≈ 0.37918288298942798287
+    @test icc(model, :subj) ≈ 0.37918288
+
 end
 
 @testset "Binomial" begin

From 204d66addc24a6f7b7d669e847a0032fc4f3c311 Mon Sep 17 00:00:00 2001
From: Phillip Alday <phillip.alday@beacon.bio>
Date: Tue, 19 Apr 2022 18:27:04 -0500
Subject: [PATCH 02/14] icc on bootstrap

---
 Project.toml             |  2 +
 src/MixedModelsExtras.jl |  1 +
 src/icc.jl               | 81 ++++++++++++++++++++++++++++++++++------
 test/Project.toml        |  5 +++
 test/icc.jl              | 27 ++++++++++++++
 5 files changed, 105 insertions(+), 11 deletions(-)

diff --git a/Project.toml b/Project.toml
index 67cb9f9..8954a90 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,11 +9,13 @@ MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 MixedModels = "3, 4"
 StatsBase = "0.33"
 StatsModels = "0.6.28"
+Tables = "1"
 julia = "1.6"
 
 [extras]
diff --git a/src/MixedModelsExtras.jl b/src/MixedModelsExtras.jl
index a42ccdd..c4bdf01 100644
--- a/src/MixedModelsExtras.jl
+++ b/src/MixedModelsExtras.jl
@@ -9,6 +9,7 @@ using MixedModels
 using Statistics
 using StatsBase
 using StatsModels
+using Tables
 
 include("icc.jl")
 include("r2.jl")
diff --git a/src/icc.jl b/src/icc.jl
index 4f0e545..cceefcd 100644
--- a/src/icc.jl
+++ b/src/icc.jl
@@ -4,17 +4,44 @@ function _group_var(vc::VarCorr, group::Symbol)
     return sum(abs2, getproperty(getproperty(vc.σρ, group), :σ))
 end
 
-function _group_var(vc::VarCorr, groups::SymbolCollection)
-    return sum(_group_var.(Ref(vc), groups))
-end
-
 function _group_var(vc::VarCorr)
     groups = collect(propertynames(vc.σρ))
     return _group_var(vc, groups)
 end
 
+# TODO rework this so that _split_by_iter is only called once
+function _group_var(vc, groups::SymbolCollection)
+    return sum(_group_var.(Ref(vc), groups))
+end
+
+# for MixedModelBootstrap
+# σtbl is boot.σs
+function _group_var(tbl, group::Symbol)
+    vals = map(values(_split_by_iter(tbl))) do iter
+        return sum(Tables.rows(iter)) do row
+            return row.group == group ? abs2(row.σ) : 0
+        end
+    end
+    return vals
+end
+
+function _group_var(tbl)
+    groups = unique((row.group for row in Tables.rows(tbl)))
+    return _group_var(tbl, groups)
+end
+
+function _split_by_iter(tbl)
+    d = Dict{Int,Vector}()
+    for row in Tables.rows(tbl)
+        vv = get!(d, row.iter, Vector{Any}())
+        push!(vv, row)
+    end
+    return d
+end
+
 """
     icc(model::MixedModel, [groups])
+    icc(boot::MixedModelBootstrap, [family], [groups])
 
 Compute the intra-class correlation coefficient (ICC) for a mixed model.
 
@@ -27,6 +54,12 @@ A single `group` can be specified as a Symbol, e.g. `:subj` or a number of group
 can be specified as an array: `[:subj, :item]`. If no `groups` are specified, then
 all grouping variables are used.
 
+When a `MixedModelBootstrap` is passed, a vector of ICC values for each
+bootstrap iteration is returned. Because `MixedModelBootstrap` does not
+store the associated model family for generalized linear mixed models,
+the family must be specified (e.g., `Bernoulli()`, `Poisson()`). A shortest
+coverage interval can be computed with `MixedModels.shortestcovint`.
+
 !!! note
     The value returned here is sometimes called the "adjusted ICC" and does not take
     the variance of the fixed effects into account (the "conditional ICC").
@@ -40,17 +73,15 @@ function icc(model::GeneralizedLinearMixedModel,
     dispersion_parameter(model) &&
         throw(ArgumentError("GLMMs with dispersion parameters are not currently supported."))
 
-    if model.resp.d isa Union{Binomial,Bernoulli}
-        σ²res = π^2 / 3
-    elseif model.resp.d isa Poisson
-        σ²res = 1.0
-    else
-        throw(ArgumentError("Family $(typeof(model.resp.d)) currently unsupported, please file an issue."))
-    end
+    σ²res = _residual_variance(model.resp.d)
 
     return _icc(VarCorr(model), groups, σ²res)
 end
 
+_residual_variance(::Union{Binomial,Bernoulli}) = π^2 / 3
+_residual_variance(::Poisson) = 1.0
+_residual_variance(::Any) = throw(ArgumentError("Family $(typeof(d)) currently unsupported, please file an issue."))
+
 function icc(model::LinearMixedModel,
              groups::Union{Symbol,SymbolCollection})
     return _icc(VarCorr(model), groups, varest(model))
@@ -60,6 +91,34 @@ icc(model::MixedModel) = icc(model, fnames(model))
 
 function _icc(vc::VarCorr, groups::Union{Symbol,SymbolCollection}, σ²res)
     σ²_α = _group_var(vc, groups) # random effect(s) in numerator
+    # TODO: don't recompute this for groups already in the above
     σ² = σ²res + _group_var(vc)
     return σ²_α / σ²
 end
+
+# TODO: upstream
+# fnames(boot::MixedModelBootstrap) = propertynames(boot.fcnames)
+icc(boot::MixedModelBootstrap) = icc(boot, propertynames(boot.fcnames))
+icc(boot::MixedModelBootstrap, family) = icc(boot, family, propertynames(boot.fcnames))
+
+function icc(boot::MixedModelBootstrap,
+            groups::Union{Symbol,SymbolCollection})
+    any(ismissing, boot.σ) &&
+        throw(ArgumentError("Bootstrapping GLMM requires specifying the family."))
+    return _icc(boot.σs, groups, abs2.(boot.σ))
+end
+
+
+function icc(boot::MixedModelBootstrap, family,
+             groups::Union{Symbol,SymbolCollection})
+    σ²res = _residual_variance(family)
+    return _icc(boot.σs, groups, σ²res)
+end
+
+
+function _icc(tbl, groups::Union{Symbol,SymbolCollection}, σ²res)
+    σ²_α = _group_var(tbl, groups) # random effect(s) in numerator
+    # TODO: don't recompute this for groups already in the above
+    σ² = σ²res .+ _group_var(tbl)
+    return σ²_α ./ σ²
+end
diff --git a/test/Project.toml b/test/Project.toml
index 6fe61e0..4f35656 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,11 @@ GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+DataFrames = "1"
+StableRNGs = "1"
\ No newline at end of file
diff --git a/test/icc.jl b/test/icc.jl
index 9bccacf..002d4c6 100644
--- a/test/icc.jl
+++ b/test/icc.jl
@@ -2,6 +2,7 @@ using DataFrames
 using MixedModels
 using MixedModels: dataset
 using MixedModelsExtras
+using StableRNGs
 using Statistics
 using Test
 
@@ -13,6 +14,23 @@ progress = false
     @test icc(model, :subj) == icc(model, [:subj]) == icc(model)
     @test icc(model, :subj) ≈ 0.37918288
 
+    formula = @formula(rt_trunc ~ 1 + spkr * prec * load +
+                                (1 + spkr | subj) +
+                                (1 | item))
+    model = fit(MixedModel, formula, dataset(:kb07); progress)
+    @test icc(model, :subj) + icc(model, :item) ≈ icc(model)
+
+    @testset "bootstrap" begin
+        boot = parametricbootstrap(StableRNG(42), 100, model; hide_progress=!progress)
+        iccboot_subj = icc(boot, :subj)
+        iccboot_item = icc(boot, :item)
+        @test iccboot_subj + iccboot_item ≈ icc(boot)
+
+        ci_subj = shortestcovint(iccboot_subj)
+        ci_item = shortestcovint(iccboot_item)
+        @test first(ci_subj) < icc(model, :subj) < last(ci_subj)
+        @test first(ci_item) < icc(model, :item) < last(ci_item)
+    end
 end
 
 @testset "Binomial" begin
@@ -32,6 +50,15 @@ end
                    contra, Binomial(); fast=true, wts=ones(length(contra.dist)), progress)
     # Bernoullis are a special case of binomial, so make sure they give the same answer
     @test icc(modelbern, Symbol("urban & dist")) ≈ icc(modelbin, Symbol("urban & dist"))
+
+    @testset "bootstrap" begin
+        boot = parametricbootstrap(StableRNG(42), 100, modelbern; hide_progress=!progress)
+        @test_throws ArgumentError icc(boot)
+        iccboot = icc(boot, Bernoulli())
+        ci = shortestcovint(iccboot)
+        @test first(ci) < icc(modelbern) < last(ci)
+        @test iccboot ≈ icc(boot, Bernoulli(), Symbol("urban & dist"))
+    end
 end
 
 @testset "Poisson" begin

From 276f995014f863cf04f628a8086476be3e08a488 Mon Sep 17 00:00:00 2001
From: Phillip Alday <phillip.alday@beacon.bio>
Date: Tue, 19 Apr 2022 18:28:19 -0500
Subject: [PATCH 03/14] patch bump

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 8954a90..691ae5a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MixedModelsExtras"
 uuid = "781a26e1-49f4-409a-8f4c-c3159d78c17e"
 authors = ["Phillip Alday <me@phillipalday.com> and contributors"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From b7cfaa1809916561445c3df0bee3de8a7fb49ec1 Mon Sep 17 00:00:00 2001
From: Phillip Alday <palday@users.noreply.github.com>
Date: Tue, 19 Apr 2022 18:30:59 -0500
Subject: [PATCH 04/14] Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 src/icc.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/icc.jl b/src/icc.jl
index cceefcd..56c290b 100644
--- a/src/icc.jl
+++ b/src/icc.jl
@@ -80,7 +80,9 @@ end
 
 _residual_variance(::Union{Binomial,Bernoulli}) = π^2 / 3
 _residual_variance(::Poisson) = 1.0
-_residual_variance(::Any) = throw(ArgumentError("Family $(typeof(d)) currently unsupported, please file an issue."))
+function _residual_variance(::Any)
+    throw(ArgumentError("Family $(typeof(d)) currently unsupported, please file an issue."))
+end
 
 function icc(model::LinearMixedModel,
              groups::Union{Symbol,SymbolCollection})
@@ -102,20 +104,18 @@ icc(boot::MixedModelBootstrap) = icc(boot, propertynames(boot.fcnames))
 icc(boot::MixedModelBootstrap, family) = icc(boot, family, propertynames(boot.fcnames))
 
 function icc(boot::MixedModelBootstrap,
-            groups::Union{Symbol,SymbolCollection})
+             groups::Union{Symbol,SymbolCollection})
     any(ismissing, boot.σ) &&
         throw(ArgumentError("Bootstrapping GLMM requires specifying the family."))
     return _icc(boot.σs, groups, abs2.(boot.σ))
 end
 
-
 function icc(boot::MixedModelBootstrap, family,
              groups::Union{Symbol,SymbolCollection})
     σ²res = _residual_variance(family)
     return _icc(boot.σs, groups, σ²res)
 end
 
-
 function _icc(tbl, groups::Union{Symbol,SymbolCollection}, σ²res)
     σ²_α = _group_var(tbl, groups) # random effect(s) in numerator
     # TODO: don't recompute this for groups already in the above

From f08fddd4915024664fe723814d44a6a21a7e7c6b Mon Sep 17 00:00:00 2001
From: Phillip Alday <palday@users.noreply.github.com>
Date: Tue, 19 Apr 2022 18:31:34 -0500
Subject: [PATCH 05/14] Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 test/icc.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/icc.jl b/test/icc.jl
index 002d4c6..57bf9de 100644
--- a/test/icc.jl
+++ b/test/icc.jl
@@ -15,8 +15,8 @@ progress = false
     @test icc(model, :subj) ≈ 0.37918288
 
     formula = @formula(rt_trunc ~ 1 + spkr * prec * load +
-                                (1 + spkr | subj) +
-                                (1 | item))
+                                  (1 + spkr | subj) +
+                                  (1 | item))
     model = fit(MixedModel, formula, dataset(:kb07); progress)
     @test icc(model, :subj) + icc(model, :item) ≈ icc(model)
 

From a7dce025bd5d61941db49684eb6c1f403cc52a0d Mon Sep 17 00:00:00 2001
From: Phillip Alday <phillip.alday@beacon.bio>
Date: Tue, 19 Apr 2022 18:32:02 -0500
Subject: [PATCH 06/14] newline

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 4f35656..7ef2835 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -11,4 +11,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 DataFrames = "1"
-StableRNGs = "1"
\ No newline at end of file
+StableRNGs = "1"

From 3ffe49ba8b25a9fd082d3297bb5b2d28086f49cc Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 11 Sep 2023 16:11:52 -0500
Subject: [PATCH 07/14] bootstrapped LRT

---
 Project.toml             |  7 +++++--
 src/MixedModelsExtras.jl |  6 ++++++
 src/bootstrap.jl         | 37 +++++++++++++++++++++++++++++++++++++
 test/bootstrap.jl        | 10 ++++++++++
 4 files changed, 58 insertions(+), 2 deletions(-)
 create mode 100644 src/bootstrap.jl
 create mode 100644 test/bootstrap.jl

diff --git a/Project.toml b/Project.toml
index 777dee9..409cb19 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,15 +6,18 @@ version = "1.1.0"
 [deps]
 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"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-MixedModels = "3, 4"
+MixedModels = "4"
+ProgressMeter = "1"
 StatsBase = "0.33, 0.34"
-StatsModels = "0.6.28, 0.7"
+StatsModels = "0.7"
 Tables = "1"
 julia = "1.6"
 
diff --git a/src/MixedModelsExtras.jl b/src/MixedModelsExtras.jl
index a66dc21..a73a86b 100644
--- a/src/MixedModelsExtras.jl
+++ b/src/MixedModelsExtras.jl
@@ -2,11 +2,14 @@ module MixedModelsExtras
 
 using LinearAlgebra
 using MixedModels
+using Random
 using Statistics
 using StatsBase
 using StatsModels
 using Tables
 
+using MixedModels: replicate
+
 include("icc.jl")
 export icc
 
@@ -22,4 +25,7 @@ export partial_fitted
 include("shrinkage.jl")
 export shrinkagenorm, shrinkagetables
 
+include("bootstrap.jl")
+export bootstrap_lrtest
+
 end
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
new file mode 100644
index 0000000..49635b7
--- /dev/null
+++ b/src/bootstrap.jl
@@ -0,0 +1,37 @@
+
+function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, m1::MixedModel; optsum_overrides=(;), progress=true)
+    models = [m0, m1]
+    dofs = dof.(models)
+    formulas = string.(formula.(models))
+    ord = sortperm(dofs)
+    dofs = dofs[ord]
+    formulas = formulas[ord]
+    devs = deviance.(models)[ord]
+    dofdiffs = diff(dofs)
+    devdiffs = .-(diff(devs))
+  
+    m0 = deepcopy(m0)
+    m1 = deepcopy(m1)
+    for (key, val) in pairs(optsum_overrides)
+        setfield!(m0.optsum, key, val)
+        setfield!(m1.optsum, key, val)
+    end
+    nulldist = replicate(n; hide_progress=!progress) do
+        simulate!(rng, m0)
+        refit!(m0; progress=false)
+        refit!(m1, response(m0))
+        return deviance(m1) - deviance(m0)
+    end
+    pvals = map(devdiffs) do dev
+        if dev > 0
+            mean(>(dev), nulldist)
+        else
+            NaN
+        end
+    end
+
+    return MixedModels.LikelihoodRatioTest(formulas,
+                                           (dof=dofs, deviance=devs),
+                                           (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),
+                                           first(models) isa LinearMixedModel)
+end
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
new file mode 100644
index 0000000..14194c9
--- /dev/null
+++ b/test/bootstrap.jl
@@ -0,0 +1,10 @@
+using MixedModels
+using MixedModelsExtras
+using Random
+using Test
+
+using MixedModels: likelihoodratiotest
+
+sleepstudy = MixedModels.dataset(:sleepstudy)
+fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), sleepstudy)
+fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)

From e72d16f5c92001fb07f6cadf2c547ed4cd9b5c25 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Mon, 11 Sep 2023 16:30:19 -0500
Subject: [PATCH 08/14] varargs

---
 src/bootstrap.jl  | 25 ++++++++++++++++---------
 test/bootstrap.jl |  1 +
 2 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 49635b7..3f43fa4 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -1,6 +1,9 @@
 
-function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, m1::MixedModel; optsum_overrides=(;), progress=true)
-    models = [m0, m1]
+function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...; optsum_overrides=(;), progress=true)
+    m0 = deepcopy(m0)
+    ms = [deepcopy(m) for m in ms]
+    
+    models = [m0; ms...]
     dofs = dof.(models)
     formulas = string.(formula.(models))
     ord = sortperm(dofs)
@@ -10,21 +13,25 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, m1::Mixe
     dofdiffs = diff(dofs)
     devdiffs = .-(diff(devs))
   
-    m0 = deepcopy(m0)
-    m1 = deepcopy(m1)
     for (key, val) in pairs(optsum_overrides)
         setfield!(m0.optsum, key, val)
-        setfield!(m1.optsum, key, val)
+        for m in ms
+            setfield!(m.optsum, key, val)
+        end
     end
     nulldist = replicate(n; hide_progress=!progress) do
         simulate!(rng, m0)
         refit!(m0; progress=false)
-        refit!(m1, response(m0))
-        return deviance(m1) - deviance(m0)
+        for m in ms
+            refit!(m, response(m0); progress=false)
+        end
+        return [deviance(m) for m in models]
     end
-    pvals = map(devdiffs) do dev
+    nulldist = stack(nulldist; dims=1)
+    nulldist = -1 .* diff(nulldist; dims=2)
+    pvals = map(enumerate(devdiffs)) do (idx, dev)
         if dev > 0
-            mean(>(dev), nulldist)
+            mean(>=(dev), view(nulldist, :, idx))
         else
             NaN
         end
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index 14194c9..dc42feb 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -8,3 +8,4 @@ using MixedModels: likelihoodratiotest
 sleepstudy = MixedModels.dataset(:sleepstudy)
 fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), sleepstudy)
 fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)
+fmzc = fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days|subj)), sleepstudy)

From 2d4d2ade4bb0910cfe285daa3611da454db3e5a7 Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Tue, 12 Sep 2023 02:09:00 -0500
Subject: [PATCH 09/14] tweak

---
 src/bootstrap.jl | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 3f43fa4..400009e 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -1,7 +1,11 @@
-
-function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...; optsum_overrides=(;), progress=true)
-    m0 = deepcopy(m0)
-    ms = [deepcopy(m) for m in ms]
+"""
+    bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
+                     optsum_overrides=(;), progress=true)
+"""
+function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...; 
+                          optsum_overrides=(;), progress=true)
+    y0 = response(m0)
+    ys = [response(m) for m in ms]
     
     models = [m0; ms...]
     dofs = dof.(models)
@@ -37,6 +41,15 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
         end
     end
 
+    # 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
+
     return MixedModels.LikelihoodRatioTest(formulas,
                                            (dof=dofs, deviance=devs),
                                            (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),

From 0c68d517deecf98668fce448150f8898be9b3ace Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Tue, 12 Sep 2023 02:21:53 -0500
Subject: [PATCH 10/14] docstring

---
 src/bootstrap.jl | 104 +++++++++++++++++++++++++++++------------------
 1 file changed, 64 insertions(+), 40 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 400009e..7007821 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -1,55 +1,79 @@
 """
     bootstrap_lrtest(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 deviance is differencce 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_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...; 
                           optsum_overrides=(;), progress=true)
     y0 = response(m0)
     ys = [response(m) for m in ms]
-    
-    models = [m0; ms...]
-    dofs = dof.(models)
-    formulas = string.(formula.(models))
-    ord = sortperm(dofs)
-    dofs = dofs[ord]
-    formulas = formulas[ord]
-    devs = deviance.(models)[ord]
-    dofdiffs = diff(dofs)
-    devdiffs = .-(diff(devs))
-  
-    for (key, val) in pairs(optsum_overrides)
-        setfield!(m0.optsum, key, val)
-        for m in ms
-            setfield!(m.optsum, key, val)
+    try
+        models = [m0; ms...]
+        dofs = dof.(models)
+        formulas = string.(formula.(models))
+        ord = sortperm(dofs)
+        dofs = dofs[ord]
+        formulas = formulas[ord]
+        devs = deviance.(models)[ord]
+        dofdiffs = diff(dofs)
+        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
-    end
-    nulldist = replicate(n; hide_progress=!progress) do
-        simulate!(rng, m0)
-        refit!(m0; progress=false)
-        for m in ms
-            refit!(m, response(m0); progress=false)
+        nulldist = replicate(n; hide_progress=!progress) do
+            simulate!(rng, m0)
+            refit!(m0; progress=false)
+            for m in ms
+                refit!(m, response(m0); progress=false)
+            end
+            return [deviance(m) for m in models]
         end
-        return [deviance(m) for m in models]
-    end
-    nulldist = stack(nulldist; dims=1)
-    nulldist = -1 .* diff(nulldist; dims=2)
-    pvals = map(enumerate(devdiffs)) do (idx, dev)
-        if dev > 0
-            mean(>=(dev), view(nulldist, :, idx))
-        else
-            NaN
+        nulldist = stack(nulldist; dims=1)
+        nulldist = -1 .* diff(nulldist; dims=2)
+        pvals = map(enumerate(devdiffs)) do (idx, dev)
+            if dev > 0
+                mean(>=(dev), view(nulldist, :, idx))
+            else
+                NaN
+            end
+        end
+    # catch ex
+    #     rethrow(ex)
+    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
-
-    # 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
-
     return MixedModels.LikelihoodRatioTest(formulas,
                                            (dof=dofs, deviance=devs),
                                            (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),

From 2eb6c6d74e98bcaeb0fff4d503a31e0b5b0b5dbf Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Tue, 12 Sep 2023 02:22:34 -0500
Subject: [PATCH 11/14] YAS

---
 src/bootstrap.jl  | 11 ++++++-----
 test/bootstrap.jl |  7 ++++---
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 7007821..95e2f5d 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -24,7 +24,7 @@ original models is compared against this null distribution.
 
 This functionality may be deprecated in the future in favor of `StatsModels.lrtest`.
 """
-function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...; 
+function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
                           optsum_overrides=(;), progress=true)
     y0 = response(m0)
     ys = [response(m) for m in ms]
@@ -38,7 +38,7 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
         devs = deviance.(models)[ord]
         dofdiffs = diff(dofs)
         devdiffs = .-(diff(devs))
-      
+
         for (key, val) in pairs(optsum_overrides)
             setfield!(m0.optsum, key, val)
             for m in ms
@@ -62,8 +62,8 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
                 NaN
             end
         end
-    # catch ex
-    #     rethrow(ex)
+        # catch ex
+        #     rethrow(ex)
     finally
         # restore the original fits
         if progress
@@ -76,6 +76,7 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
     end
     return MixedModels.LikelihoodRatioTest(formulas,
                                            (dof=dofs, deviance=devs),
-                                           (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),
+                                           (dofdiff=dofdiffs, deviancediff=devdiffs,
+                                            pvalues=pvals),
                                            first(models) isa LinearMixedModel)
 end
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index dc42feb..71df1a7 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -6,6 +6,7 @@ using Test
 using MixedModels: likelihoodratiotest
 
 sleepstudy = MixedModels.dataset(:sleepstudy)
-fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), sleepstudy)
-fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)
-fmzc = fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days|subj)), sleepstudy)
+fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | subj)), sleepstudy)
+fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), sleepstudy)
+fmzc = fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days | subj)),
+           sleepstudy)

From 50980580fbe1cd9de5605dec3c5826e9efc4cb1f Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Tue, 12 Sep 2023 03:02:19 -0500
Subject: [PATCH 12/14] try-catch instead of copy?

---
 src/MixedModelsExtras.jl |  2 +-
 src/bootstrap.jl         | 20 ++++++++++++++------
 test/bootstrap.jl        | 12 +++++++++---
 3 files changed, 24 insertions(+), 10 deletions(-)

diff --git a/src/MixedModelsExtras.jl b/src/MixedModelsExtras.jl
index a73a86b..8e610c2 100644
--- a/src/MixedModelsExtras.jl
+++ b/src/MixedModelsExtras.jl
@@ -26,6 +26,6 @@ include("shrinkage.jl")
 export shrinkagenorm, shrinkagetables
 
 include("bootstrap.jl")
-export bootstrap_lrtest
+export bootstrap_lrt
 
 end
diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 95e2f5d..9bcc533 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -1,6 +1,6 @@
 """
-    bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
-                     optsum_overrides=(;), progress=true)
+    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.
 
@@ -24,10 +24,18 @@ original models is compared against this null distribution.
 
 This functionality may be deprecated in the future in favor of `StatsModels.lrtest`.
 """
-function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
+function bootstrap_lrt(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
                           optsum_overrides=(;), progress=true)
     y0 = response(m0)
     ys = [response(m) for m in ms]
+    local models
+    local devs
+    local dofs
+    local dofs
+    local formulas
+    local dofdiffs
+    local devdiffs
+    local pvals
     try
         models = [m0; ms...]
         dofs = dof.(models)
@@ -54,7 +62,7 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
             return [deviance(m) for m in models]
         end
         nulldist = stack(nulldist; dims=1)
-        nulldist = -1 .* diff(nulldist; dims=2)
+        nulldist = diff(nulldist; dims=2)
         pvals = map(enumerate(devdiffs)) do (idx, dev)
             if dev > 0
                 mean(>=(dev), view(nulldist, :, idx))
@@ -62,8 +70,8 @@ function bootstrap_lrtest(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::Mixe
                 NaN
             end
         end
-        # catch ex
-        #     rethrow(ex)
+    catch ex
+        rethrow(ex)
     finally
         # restore the original fits
         if progress
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index 71df1a7..e372a76 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -5,8 +5,14 @@ using Test
 
 using MixedModels: likelihoodratiotest
 
+progress = false
 sleepstudy = MixedModels.dataset(:sleepstudy)
-fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | subj)), sleepstudy)
-fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), 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)
+           sleepstudy; progress)
+
+lrt = likelihoodratiotest(fm0, fm1, fmzc)
+boot = bootstrap_lrt(MersenneTwister(42), 1000, fm0, fm1, fmzc)

From 540285cd9ca4dec9a850bfc1fc96031efc2d363d Mon Sep 17 00:00:00 2001
From: Phillip Alday <me@phillipalday.com>
Date: Tue, 12 Sep 2023 06:09:28 -0500
Subject: [PATCH 13/14] d'oh

---
 src/bootstrap.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 9bcc533..4f9e0a3 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -26,8 +26,8 @@ This functionality may be deprecated in the future in favor of `StatsModels.lrte
 """
 function bootstrap_lrt(rng::AbstractRNG, n::Integer, m0::MixedModel, ms::MixedModel...;
                           optsum_overrides=(;), progress=true)
-    y0 = response(m0)
-    ys = [response(m) for m in ms]
+    y0 = copy(response(m0))
+    ys = [copy(response(m)) for m in ms]
     local models
     local devs
     local dofs

From 07f0ca1e48360764507e18613ebc1a9782272d88 Mon Sep 17 00:00:00 2001
From: Phillip Alday <phillip.alday@beacon.bio>
Date: Thu, 14 Sep 2023 10:12:55 +0200
Subject: [PATCH 14/14] YAS

---
 src/bootstrap.jl  | 2 +-
 test/bootstrap.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/bootstrap.jl b/src/bootstrap.jl
index 4f9e0a3..7951a06 100644
--- a/src/bootstrap.jl
+++ b/src/bootstrap.jl
@@ -25,7 +25,7 @@ original models is compared against this null distribution.
 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)
+                       optsum_overrides=(;), progress=true)
     y0 = copy(response(m0))
     ys = [copy(response(m)) for m in ms]
     local models
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
index e372a76..cade49b 100644
--- a/test/bootstrap.jl
+++ b/test/bootstrap.jl
@@ -7,9 +7,9 @@ using MixedModels: likelihoodratiotest
 
 progress = false
 sleepstudy = MixedModels.dataset(:sleepstudy)
-fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | subj)), 
+fm0 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | subj)),
           sleepstudy; progress)
-fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), 
+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)