shortestcovint for MixedModelBootstrap objects (#484)

* shortestcovint for MixedModelBootstrap objects

* news update

* minor version bump

* release notes text update

* showcase new method in docs

.github/workflows/TagBot.yml

@@ -15,26 +15,31 @@ jobs:
           ssh: ${{ secrets.DOCUMENTER_KEY }} # see https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#GitHub-Actions
           changelog: |
               ## {{ package }} {{ version }}
-              
+
               {% if previous_release %}
-              [NEWS file](https://github.com/JuliaStats/MixedModels.jl/blob/{{ version }}/NEWS.md). 
+              [NEWS file](https://github.com/JuliaStats/MixedModels.jl/blob/{{ version }}/NEWS.md).
               [Diff since {{ previous_release }}]({{ compare_url }})
               {% endif %}
-              
+
               {% if custom %}
               {{ custom }}
               {% endif %}
-              
+
+
+              *NB: Closed issues and pull requests are sorted temporally and so may
+              include backports to other versions or work in the development branch for
+              an upcoming breaking release. Please see the [NEWS file](https://github.com/JuliaStats/MixedModels.jl/blob/{{ version }}/NEWS.md)
+              for changes sorted by release.*
               {% if issues %}
               **Closed issues:**
               {% for issue in issues %}
               - {{ issue.title }} (#{{ issue.number }})
               {% endfor %}
               {% endif %}
-              
+
               {% if pulls %}
               **Merged pull requests:**
               {% for pull in pulls %}
               - {{ pull.title }} (#{{ pull.number }}) (@{{ pull.author.username }})
               {% endfor %}
-              {% endif %}            
+              {% endif %}

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v3.4.0 Release Notes
+========================
+* `shortestcovint` method for `MixedModelsBootstrap` [#484]
+
 MixedModels v3.3.0 Release Notes
 ========================
 * HTML and LaTeX `show` methods for `MixedModel`, `BlockDescription`,
@@ -145,3 +149,4 @@ Package dependencies
 [#449]: https://github.com/JuliaStats/MixedModels.jl/issues/449
 [#474]: https://github.com/JuliaStats/MixedModels.jl/issues/474
 [#480]: https://github.com/JuliaStats/MixedModels.jl/issues/480
+[#484]: https://github.com/JuliaStats/MixedModels.jl/issues/484

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 = "3.3.0"
+version = "3.4.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

docs/src/bootstrap.md

@@ -94,6 +94,11 @@ We generate these for all random and fixed effects:
 combine(groupby(df, [:type, :group, :names]), :value => shortestcovint => :interval)
 ```
 
+We can also generate this directly from the original bootstrap object:
+```@example Main
+DataFrame(shortestcovint(samp))
+```
+
 A value of zero for the standard deviation of the random effects is an example of a *singular* covariance.
 It is easy to detect the singularity in the case of a scalar random-effects term.
 However, it is not as straightforward to detect singularity in vector-valued random-effects terms.
@@ -116,7 +121,7 @@ first(df2, 10)
 the singularity can be exhibited as a standard deviation of zero or as a correlation of $\pm1$.
 
 ```@example Main
-combine(groupby(df2, [:type, :group, :names]), :value => shortestcovint => :interval)
+DataFrame(shortestcovint(samp2))
 ```
 
 A histogram of the estimated correlations from the bootstrap sample has a spike at `+1`.

src/bootstrap.jl

@@ -125,7 +125,7 @@ end
 """
     allpars(bsamp::MixedModelBootstrap)
 
-Return a tidy (row)table with the parameter estimates spread into columns
+Return a tidy (column)table with the parameter estimates spread into columns
 of `iter`, `type`, `group`, `name` and `value`.
 """
 function allpars(bsamp::MixedModelBootstrap{T}) where {T}
@@ -245,6 +245,48 @@ function shortestcovint(v, level = 0.95)
     (vv[i], vv[i + ilen - 1])
 end
 
+"""
+    shortestcovint(bsamp::MixedModelBootstrap, level = 0.95)
+
+Return the shortest interval containing `level` proportion for each parameter from [`bsamp.allpars`](@ref)
+"""
+function shortestcovint(bsamp::MixedModelBootstrap{T}, level = 0.95) where {T}
+    allpars = bsamp.allpars
+    pars = unique(zip(allpars.type, allpars.group, allpars.names))
+
+    colnms = (:type, :group, :names, :lower, :upper)
+    coltypes = Tuple{String, Union{Missing,String}, Union{Missing,String}, T, T}
+    # not specifying the full eltype (NamedTuple{colnms,coltypes}) leads to prettier printing
+    result = NamedTuple{colnms}[]
+    sizehint!(result, length(pars))
+
+
+    for (t, g, n) in pars
+        gidx = if ismissing(g)
+            ismissing.(allpars.group)
+        else
+            .!ismissing.(allpars.group) .& (allpars.group .== g)
+        end
+
+        nidx = if ismissing(n)
+            ismissing.(allpars.names)
+        else
+            .!ismissing.(allpars.names) .& (allpars.names .== n)
+        end
+
+        tidx = allpars.type .== t # no missings allowed here
+
+        idx = tidx .& gidx .& nidx
+
+        vv = view(allpars.value, idx)
+
+        lower, upper = shortestcovint(vv, level)
+        push!(result, (; type=t, group=g, names=n, lower=lower, upper=upper))
+    end
+
+    return result
+end
+
 """
     tidyβ(bsamp::MixedModelBootstrap)
 Return a tidy (row)table with the parameter estimates spread into columns

test/bootstrap.jl

@@ -104,11 +104,7 @@ end
         contra = dataset(:contra)
         gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=true)
         bs = parametricbootstrap(StableRNG(42), 100, gm0)
-        bsci = combine(groupby(DataFrame(bs.β), :coefname),
-                       :β => shortestcovint => :ci)
-        bsci.lower = first.(bsci.ci)
-        bsci.upper = last.(bsci.ci)
-        select!(bsci, Not(:ci))
+        bsci = filter!(:type => ==("β"), DataFrame(shortestcovint(bs)))
         ciwidth = 2 .* stderror(gm0)
         waldci = DataFrame(coef=fixefnames(gm0),
                            lower=fixef(gm0) .- ciwidth,