From 4218789b9020540c48afc2bcf4104676d66bf4a3 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 8 Mar 2021 16:09:32 +0000 Subject: [PATCH] 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 | 17 +++++++++----- NEWS.md | 5 ++++ Project.toml | 2 +- docs/src/bootstrap.md | 7 +++++- src/bootstrap.jl | 44 +++++++++++++++++++++++++++++++++++- test/bootstrap.jl | 6 +---- 6 files changed, 67 insertions(+), 14 deletions(-) diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 8f2bbc36d..6c27383e7 100644 --- a/.github/workflows/TagBot.yml +++ b/.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 %} diff --git a/NEWS.md b/NEWS.md index 90e560419..95ccbeef4 100644 --- a/NEWS.md +++ b/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 diff --git a/Project.toml b/Project.toml index 4b945c22f..e6707c467 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MixedModels" uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" author = ["Phillip Alday ", "Douglas Bates ", "Jose Bayoan Santiago Calderon "] -version = "3.3.0" +version = "3.4.0" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" diff --git a/docs/src/bootstrap.md b/docs/src/bootstrap.md index 4cb24b7f9..98241ae8a 100644 --- a/docs/src/bootstrap.md +++ b/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`. diff --git a/src/bootstrap.jl b/src/bootstrap.jl index 6af8c6fa1..cc3370a78 100644 --- a/src/bootstrap.jl +++ b/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 diff --git a/test/bootstrap.jl b/test/bootstrap.jl index 7d08423ff..c5063d565 100644 --- a/test/bootstrap.jl +++ b/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,