show method for bootstrap, profile; equal tail confint (#715)

* show method for bootstrap

* profile show method

* add support for equal tail confint on bootstrap

* news and version bump

* equaltail test

* test show methods

* test DataFrame conversion

* bump Aqua compat

* privateering

NEWS.md

@@ -1,3 +1,8 @@
+MixedModels v4.22.0 Release Notes
+==============================
+* Support for equal-tail confidence intervals for `MixedModelBootstrap`. [#715]
+* Basic `show` methods for `MixedModelBootstrap` and `MixedModelProfile`. [#715]
+
 MixedModels v4.21.0 Release Notes
 ==============================
 * Auto apply `Grouping()` to grouping variables that don't already have an explicit contrast. [#652]
@@ -476,3 +481,4 @@ Package dependencies
 [#703]: https://github.com/JuliaStats/MixedModels.jl/issues/703
 [#707]: https://github.com/JuliaStats/MixedModels.jl/issues/707
 [#709]: https://github.com/JuliaStats/MixedModels.jl/issues/709
+[#715]: https://github.com/JuliaStats/MixedModels.jl/issues/715

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>", "Jose Bayoan Santiago Calderon <jbs3hp@virginia.edu>"]
-version = "4.21.1"
+version = "4.22.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/bootstrap.jl

@@ -159,6 +159,22 @@ function Base.reduce(::typeof(vcat), v::AbstractVector{MixedModelBootstrap{T}})
         deepcopy(b1.fcnames))
 end
 
+function Base.show(io::IO, mime::MIME"text/plain", x::MixedModelBootstrap)
+    tbl = x.tbl
+    println(io, "MixedModelBootstrap with $(length(x)) samples")
+    out = NamedTuple[]
+    for col in Tables.columnnames(tbl)
+        col == :obj && continue
+        s = summarystats(Tables.getcolumn(tbl, col))
+        push!(out, (; parameter=col, s.min, s.q25, s.median, s.mean, s.q75, s.max))
+    end
+    tt = FlexTable(out)
+    # trim out the FlexTable header
+    str = last(split(sprint(show, mime, tt), "\n"; limit=2))
+    println(io, str)
+    return nothing
+end
+
 """
     parametricbootstrap([rng::AbstractRNG], nsamp::Integer, m::MixedModel{T}, ftype=T;
         β = coef(m), σ = m.σ, θ = m.θ, hide_progress=false, optsum_overrides=(;))
@@ -303,10 +319,14 @@ function allpars(bsamp::MixedModelFitCollection{T}) where {T}
 end
 
 """
-    confint(pr::MixedModelBootstrap; level::Real=0.95)
+    confint(pr::MixedModelBootstrap; level::Real=0.95, method=:shortest)
 
 Compute bootstrap confidence intervals for coefficients and variance components, with confidence level level (by default 95%).
 
+The keyword argument `method` determines whether the `:shortest`, i.e. highest density, interval is used 
+or the `:equaltail`, i.e. quantile-based, interval is used. For historical reasons, the default is `:shortest`, 
+but `:equaltail` gives the interval that is most comparable to the profile and Wald confidence intervals.
+
 !!! note
     The API guarantee is for a Tables.jl compatible table. The exact return type is an
     implementation detail and may change in a future minor release without being considered
@@ -319,7 +339,11 @@ Compute bootstrap confidence intervals for coefficients and variance components,
 
 See also [`shortestcovint`](@ref).
 """
-function StatsBase.confint(bsamp::MixedModelBootstrap{T}; level::Real=0.95) where {T}
+function StatsBase.confint(
+    bsamp::MixedModelBootstrap{T}; level::Real=0.95, method=:shortest
+) where {T}
+    method in [:shortest, :equaltail] ||
+        throw(ArgumentError("`method` must be either :shortest or :equaltail."))
     cutoff = sqrt(quantile(Chisq(1), level))
     # Creating the table is somewhat wasteful because columns are created then immediately skipped.
     tbl = Table(bsamp.tbl)
@@ -333,8 +357,13 @@ function StatsBase.confint(bsamp::MixedModelBootstrap{T}; level::Real=0.95) wher
             ),
         ),
     )
+    tails = [(1 - level) / 2, (1 + level) / 2]
     for p in par
-        l, u = shortestcovint(sort!(copyto!(v, getproperty(tbl, p))), level)
+        if method === :shortest
+            l, u = shortestcovint(sort!(copyto!(v, getproperty(tbl, p))), level)
+        else
+            l, u = quantile(getproperty(tbl, p), tails)
+        end
         push!(lower, l)
         push!(upper, u)
     end

src/profile/profile.jl

@@ -80,3 +80,11 @@ function StatsAPI.confint(pr::MixedModelProfile; level::Real=0.95)
         upper=[rev[s](cutoff) for s in syms],
     )
 end
+
+function Base.show(io::IO, mime::MIME"text/plain", pr::MixedModelProfile)
+    print(io, "MixedModelProfile -- ")
+    show(io, mime, pr.tbl)
+    return nothing
+end
+
+Tables.columns(pr::MixedModelProfile) = Tables.columns(pr.tbl)

test/Project.toml

@@ -17,6 +17,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]
-Aqua = "0.5, 0.6"
+Aqua = "0.7"
 StableRNGs = "0.1, 1"
 Suppressor = "0.2"

test/bootstrap.jl

@@ -221,14 +221,28 @@ end
     end
 end
 
-@testset "CI method comparison" begin
+@testset "show and summary" begin
     fmzc = models(:sleepstudy)[2]
-    level = 0.68
     pb = parametricbootstrap(MersenneTwister(42), 500, fmzc; hide_progress=true)
     pr = profile(fmzc)
-    ci_boot = confint(pb; level)
-    ci_wald = confint(fmzc; level)
-    ci_prof = confint(pr; level)
-    @test first(ci_boot.lower, 2) ≈ first(ci_prof.lower, 2) atol=0.5
-    @test first(ci_prof.lower, 2) ≈ first(ci_wald.lower, 2) atol=0.1
+    @test startswith(sprint(show, MIME("text/plain"), pr),
+                     "MixedModelProfile -- Table with 9 columns and 151 rows:")
+    @test startswith(sprint(show, MIME("text/plain"), pb),
+                     "MixedModelBootstrap with 500 samples\n     parameter  min        q25       median    mean      q75       max\n  ")
+
+    df = DataFrame(pr)
+    @test nrow(df) == 151
+    @test propertynames(df) == collect(propertynames(pr.tbl))
+
+    @testset "CI method comparison" begin    
+        level = 0.68
+        ci_boot_equaltail = confint(pb; level, method=:equaltail)
+        ci_boot_shortest = confint(pb; level, method=:shortest)
+        @test_throws ArgumentError confint(pb; level, method=:other)
+        ci_wald = confint(fmzc; level)
+        ci_prof = confint(pr; level)
+        @test first(ci_boot_shortest.lower, 2) ≈ first(ci_prof.lower, 2) atol=0.5
+        @test first(ci_boot_equaltail.lower, 2) ≈ first(ci_prof.lower, 2) atol=0.5
+        @test first(ci_prof.lower, 2) ≈ first(ci_wald.lower, 2) atol=0.1
+    end
 end

test/runtests.jl

@@ -1,4 +1,5 @@
 using Aqua
+using GLM
 using MixedModels
 using Test
 
@@ -12,7 +13,9 @@ import LinearAlgebra: BLAS
 @testset "Aqua" begin
     # we can't check for unbound type parameters
     # because we actually need one at one point for _same_family()
-    Aqua.test_all(MixedModels; ambiguities=false, unbound_args=false)
+    Aqua.test_all(MixedModels; ambiguities=false, unbound_args=false,
+                  # XXX TODO: upstream this piracy
+                  piracy=(;treat_as_own=[GLM.wrkresp!, Base.:|]))
 end
 
 include("utilities.jl")