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 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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