fix display and confint of bootstrap for models without dispersion parameter

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v5.0.2 Release Notes
+==============================
+- The default display and `confint` methods for bootstrap results from models without dispersion parameters has been fixed. [#861]
+
 MixedModels v5.0.1 Release Notes
 ==============================
 - Fixes a method error with `Grouping()` contrasts when used with recent CategoricalArray releases. [#860]
@@ -690,3 +694,4 @@ Package dependencies
 [#857]: https://github.com/JuliaStats/MixedModels.jl/issues/857
 [#858]: https://github.com/JuliaStats/MixedModels.jl/issues/858
 [#860]: https://github.com/JuliaStats/MixedModels.jl/issues/860
+[#861]: https://github.com/JuliaStats/MixedModels.jl/issues/861

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>"]
-version = "5.0.1"
+version = "5.0.2"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

docs/src/formula_syntax.md

@@ -41,7 +41,7 @@ If you require specific functionality from RegressionFormulae.jl, it is best to
   - `a / b` is read as "`b` is nested within `a`".
   - `a / b` expands to `a + fulldummy(a) & b`.
 - It is generally not necessary to specify nesting in the blocking variables, when the inner levels are unique across outer levels. In other words, in a study with children (`C1`, `C2`, etc. ) nested within schools (`S1`, `S2`, etc.),
-  - it is not necessary to specify the nesting when `C1` identifies a unique child across schools. In other words, intercept-only random effects terms can be written as `(1|C) + `(1|S)`.
+  - it is not necessary to specify the nesting when `C1` identifies a unique child across schools. In other words, intercept-only random effects terms can be written as `(1|C) + (1|S)`.
   - it is necessary to specify the nesting when chid identifiers are re-used across schools, e.g. `C1` refers to a child in `S1` and a different child in `S2`. In this case, the nested syntax `(1|S/C)` expands to `(1|S) + (1|S&C)`. The interaction term in the second blocking variable generates unique labels for each child across schools.
 
 

src/bootstrap.jl

@@ -355,14 +355,16 @@ function StatsBase.confint(
     tbl = Table(bsamp.tbl)
     lower = T[]
     upper = T[]
-    v = similar(tbl.σ)
-    par = sort!(
-        collect(
-            filter(
-                k -> !(startswith(string(k), 'θ') || string(k) == "obj"), propertynames(tbl)
-            ),
-        ),
-    )
+    v = similar(tbl.σ, T)
+    par = filter(collect(propertynames(tbl))) do k
+        k = string(k)
+        # σ is missing in models without a dispersion parameter
+        if k == "σ" && Missing <: eltype(tbl.σ)
+            return false
+        end
+        return !startswith(k, 'θ') && k != "obj"
+    end
+    sort!(par)
     tails = [(1 - level) / 2, (1 + level) / 2]
     for p in par
         if method === :shortest
@@ -630,6 +632,8 @@ push! `σ` times the row lengths (σs) and the inner products of normalized rows
 """
 function σρ!(v::AbstractVector{<:Union{T,Missing}}, t::LowerTriangular, σ) where {T}
     dat = t.data
+    # for models without a dispersion parameter, σ is missing, but for the math below we can treat it as 1
+    σ = coalesce(σ, one(T))
     for i in axes(dat, 1)
         ssqr = zero(T)
         for j in 1:i

test/bootstrap.jl

@@ -219,15 +219,13 @@ end
         bs = parametricbootstrap(StableRNG(42), 100, gm0; progress=false)
         # make sure we're not copying
         @test length(bs.lowerbd) == length(gm0.θ)
-        bsci = filter!(:type => ==("β"), DataFrame(shortestcovint(bs)))
-        ciwidth = 2 .* stderror(gm0)
-        waldci = DataFrame(; coef=fixefnames(gm0),
-            lower=fixef(gm0) .- ciwidth,
-            upper=fixef(gm0) .+ ciwidth)
+        bsci = confint(bs)
+        waldci = confint(gm0)
 
         # coarse tolerances because we're not doing many bootstrap samples
-        @test all(isapprox.(bsci.lower, waldci.lower; atol=0.5))
-        @test all(isapprox.(bsci.upper, waldci.upper; atol=0.5))
+        # end-1 because the bootstrap CIs include the variance component
+        @test all(isapprox.(collect(bsci.lower)[1:end-1], collect(waldci.lower); atol=0.5))
+        @test all(isapprox.(collect(bsci.upper)[1:end-1], collect(waldci.upper); atol=0.5))
 
         σbar = mean(MixedModels.tidyσs(bs)) do x
             x.σ