diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 27ab1c5..d5f9834 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,12 +22,12 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - version: [1.8] + version: [1.9] arch: [x64] os: [ubuntu-22.04] # macos-10.15, windows-2019 steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v4 - name: Julia Setup uses: julia-actions/setup-julia@v1 with: @@ -44,7 +44,7 @@ jobs: env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - name: Percy Upload - if: ${{ matrix.version == '1.8' }} + if: ${{ matrix.version == '1.9' }} run: | ls ./test/output/ # useful for debugging npx @percy/cli upload ./test/output @@ -53,4 +53,6 @@ jobs: - name: Coverage Process uses: julia-actions/julia-processcoverage@v1 - name: Coverage Upload - uses: codecov/codecov-action@v3 + uses: codecov/codecov-action@v4 + with: + token: ${{ secrets.CODECOV_TOKEN }} diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml index 50f7b3d..4a0b00e 100644 --- a/.github/workflows/documenter.yml +++ b/.github/workflows/documenter.yml @@ -25,7 +25,7 @@ jobs: # Run on push's or non-draft PRs if: (github.event_name == 'push') || (github.event.pull_request.draft == false) steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: 1 diff --git a/.github/workflows/style.yml b/.github/workflows/style.yml index 42b467f..2ebd63c 100644 --- a/.github/workflows/style.yml +++ b/.github/workflows/style.yml @@ -18,7 +18,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - julia-version: [1.8] + julia-version: [1.9] steps: - uses: julia-actions/setup-julia@latest with: diff --git a/Project.toml b/Project.toml index 5bd903b..0e08989 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MixedModelsMakie" uuid = "b12ae82c-6730-437f-aff9-d2c38332a376" authors = ["Phillip Alday ", "Douglas Bates ", "contributors"] -version = "0.3.28" +version = "0.4.0" [deps] BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" @@ -19,11 +19,11 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" [compat] BSplineKit = "0.15, 0.16, 0.17" DataFrames = "1" -Distributions = "0.21, 0.22, 0.23, 0.24, 0.25" +Distributions = "0.25" KernelDensity = "0.6.3" Makie = "0.20" MixedModels = "4.14" PrecompileTools = "1" SpecialFunctions = "1, 2" -StatsBase = "0.31, 0.32, 0.33, 0.34" -julia = "1.8" +StatsBase = "0.33, 0.34" +julia = "1.9" diff --git a/src/MixedModelsMakie.jl b/src/MixedModelsMakie.jl index e1e1b50..870e42c 100644 --- a/src/MixedModelsMakie.jl +++ b/src/MixedModelsMakie.jl @@ -36,6 +36,10 @@ export RanefInfo, zetaplot, zetaplot! +# from https://github.com/MakieOrg/Makie.jl/issues/2992 +const Indexable = Union{Makie.Figure,Makie.GridLayout,Makie.GridPosition, + Makie.GridSubposition} + include("utilities.jl") include("shrinkage.jl") include("caterpillar.jl") @@ -48,7 +52,7 @@ include("recipes.jl") @setup_workload begin model = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days | subj)), - MixedModels.dataset(:sleepstudy)) + MixedModels.dataset(:sleepstudy); progress=false) @compile_workload begin caterpillar(model) coefplot(model) diff --git a/src/coefplot.jl b/src/coefplot.jl index 59a653c..60a02be 100644 --- a/src/coefplot.jl +++ b/src/coefplot.jl @@ -1,65 +1,59 @@ """ - coefplot(x::MixedModel; - conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...) - coefplot(x::MixedModelBootstrap; - conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...) + coefplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure + coefplot!(fig::$(Indexable), x::Union{MixedModel,MixedModelBootstrap}; + kwargs...) + coefplot!(ax::Axis, Union{MixedModel,MixedModelBootstrap}; + conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...) Create a coefficient plot of the fixed-effects and associated confidence intervals. -!!! note - This functionality is implemented using [Makie recipes](https://makie.juliaplots.org/v0.15.0/recipes.html) - and thus there are also additional auto-generated methods for `coefplot` and `coefplot!` that may be useful - when constructing more complex figures. +Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded. + +`attributes` are passed onto `scatter!` and `errorbars!`. + +The mutating methods return the original object. + +!!! note + Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) + are excluded. """ -function coefplot(x::Union{MixedModel,MixedModelBootstrap}; - conf_level=0.95, - vline_at_zero=true, - show_intercept=true, - attributes...) +function coefplot(x::Union{MixedModel,MixedModelBootstrap}; show_intercept=true, kwargs...) # need to guarantee a min height of 150 fig = Figure(; size=(640, max(150, 75 * _npreds(x; show_intercept)))) - ax = Axis(fig[1, 1]) - pl = coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, attributes...) - return Makie.FigureAxisPlot(fig, ax, pl) + coefplot!(fig, x; kwargs...) + return fig end -@recipe(CoefPlot, x) do scene - return Attributes(; conf_level=0.95, vline_at_zero=true, show_intercept=true) +"""$(@doc coefplot)""" +function coefplot!(fig::Indexable, x::Union{MixedModel,MixedModelBootstrap}; kwargs...) + ax = Axis(fig[1, 1]) + coefplot!(ax, x; kwargs...) + return fig end -function Makie.plot!(ax::Axis, P::Type{<:CoefPlot}, allattrs::Makie.Attributes, x) - plot = Makie.plot!(ax.scene, P, allattrs, x) +"""$(@doc coefplot)""" +function coefplot!(ax::Axis, x::Union{MixedModel,MixedModelBootstrap}; + conf_level=0.95, + vline_at_zero=true, + show_intercept=true, + attributes...) + ci = confint_table(x, conf_level; show_intercept) + y = nrow(ci):-1:1 + xvals = ci.estimate + xlabel = @sprintf "Estimate and %g%% confidence interval" conf_level * 100 + + attributes = merge((; xlabel), attributes) + ax.xlabel = attributes.xlabel + + scatter!(ax, xvals, y; attributes...) + errorbars!(ax, xvals, y, xvals .- ci.lower, ci.upper .- xvals; + direction=:x, attributes...) + vline_at_zero && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) - if haskey(allattrs, :title) - ax.title = allattrs.title[] - end - if haskey(allattrs, :xlabel) - ax.xlabel = allattrs.xlabel[] - else - ax.xlabel = @sprintf "Estimate and %g%% confidence interval" (allattrs.conf_level[] * - 100) - end - if haskey(allattrs, :ylabel) - ax.ylabel = allattrs.ylabel[] - end reset_limits!(ax) - show_intercept = allattrs.show_intercept[] cn = _coefnames(x; show_intercept) nticks = _npreds(x; show_intercept) ax.yticks = (nticks:-1:1, cn) ylims!(ax, 0, nticks + 1) - allattrs.vline_at_zero[] && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) - return plot -end - -function Makie.plot!(plot::CoefPlot{<:Tuple{Union{MixedModel,MixedModelBootstrap}}}) - model_or_boot, conf_level = plot[1][], plot.conf_level[] - ci = confint_table(model_or_boot, conf_level) - plot.show_intercept[] || filter!(:coefname => !=("(Intercept)"), ci) - y = nrow(ci):-1:1 - xvals = ci.estimate - scatter!(plot, xvals, y) - errorbars!(plot, xvals, y, xvals .- ci.lower, ci.upper .- xvals; direction=:x) - - return plot + return ax end diff --git a/src/ridge.jl b/src/ridge.jl index 40fc70d..d2ad163 100644 --- a/src/ridge.jl +++ b/src/ridge.jl @@ -1,15 +1,9 @@ - -@recipe(RidgePlot, x) do scene - return Attributes(; - conf_level=0.95, - vline_at_zero=true, - show_intercept=true) -end - """ - ridgeplot(x::MixedModelBootstrap; - conf_level=0.95, vline_at_zero=true, show_intercept=true, - attributes...) + ridgeplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure + ridgeplot!(fig::$(Indexable), x::Union{MixedModel,MixedModelBootstrap}; + kwargs...) + ridgeplot!(ax::Axis, Union{MixedModel,MixedModelBootstrap}; + conf_level=0.95, vline_at_zero=true, show_intercept=true, attributes...) Create a ridge plot for the bootstrap samples of the fixed effects. @@ -18,84 +12,89 @@ Densities are normalized so that the maximum density is always 1. The highest density interval corresponding to `conf_level` is marked with a bar at the bottom of each density. Setting `conf_level=missing` removes the markings for the highest density interval. -!!! note - This functionality is implemented using [Makie recipes](https://makie.juliaplots.org/v0.15.0/recipes.html) - and thus there are also additional auto-generated methods for `ridgeplot` and `ridgeplot!` that may be useful - when constructing more complex figures. +`attributes` are passed onto [`coefplot`](@ref), `band!` and `lines!`. + +The mutating methods return the original object. + +!!! note + Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) + are excluded. """ -function ridgeplot(x::MixedModelBootstrap; - conf_level=0.95, - vline_at_zero=true, - show_intercept=true, - attributes...) +function ridgeplot(x::MixedModelBootstrap; show_intercept=true, kwargs...) # need to guarantee a min height of 200 fig = Figure(; size=(640, max(200, 100 * _npreds(x; show_intercept)))) - ax = Axis(fig[1, 1]) - if !ismissing(conf_level) - pl = coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, color=:black, - attributes...) - end - pl = ridgeplot!(ax, x; vline_at_zero, conf_level, show_intercept, attributes...) - return Makie.FigureAxisPlot(fig, ax, pl) + return ridgeplot!(fig, x; show_intercept, kwargs...) end -function Makie.plot!(ax::Axis, P::Type{<:RidgePlot}, allattrs::Makie.Attributes, x) - plot = Makie.plot!(ax.scene, P, allattrs, x) +"""$(@doc ridgeplot)""" +function ridgeplot!(fig::Indexable, x::MixedModelBootstrap; kwargs...) + ax = Axis(fig[1, 1]) + ridgeplot!(ax, x; kwargs...) + return fig +end - if haskey(allattrs, :title) - ax.title = allattrs.title[] - end - if haskey(allattrs, :xlabel) - ax.xlabel = allattrs.xlabel[] +""" + _color(s::Symbol) + _color(p::Pair) + +Extract the color part out of either a color name or a `(color, alpha)` pair. +""" +_color(s) = s +_color(p::Pair) = first(p) + +"""$(@doc ridgeplot)""" +function ridgeplot!(ax::Axis, x::MixedModelBootstrap; + conf_level=0.95, + vline_at_zero=true, + show_intercept=true, + attributes...) + xlabel = if !ismissing(conf_level) + @sprintf "Normalized bootstrap density and %g%% confidence interval" (conf_level * + 100) else - lab = if !ismissing(allattrs.conf_level[]) - @sprintf "Normalized bootstrap density and %g%% confidence interval" (allattrs.conf_level[] * - 100) - else - "Normalized bootstrap density" - end - ax.xlabel = lab - end - if haskey(allattrs, :ylabel) - ax.ylabel = allattrs.ylabel[] + "Normalized bootstrap density" end - reset_limits!(ax) - show_intercept = allattrs.show_intercept[] - # check conf_level so that we don't double print - # if coefplot took care of it for us - if ismissing(allattrs.conf_level[]) - cn = _coefnames(x; show_intercept) - nticks = _npreds(x; show_intercept) - ax.yticks = (nticks:-1:1, cn) - ylims!(ax, 0, nticks + 1) - allattrs.vline_at_zero[] && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) + if !ismissing(conf_level) + coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, color=:black, + attributes...) end - return plot -end -function Makie.plot!(plot::RidgePlot{<:Tuple{MixedModelBootstrap}}) - boot, conf_level = plot[1][], plot.conf_level[] - df = transform!(DataFrame(boot.β), :coefname => ByRow(string) => :coefname) - plot.show_intercept[] || filter!(:coefname => !=("(Intercept)"), df) + attributes = merge((; xlabel, color=:black), attributes) + band_attributes = merge(attributes, (; color=(_color(attributes.color), 0.3))) + + ax.xlabel = attributes.xlabel + + df = transform!(DataFrame(x.β), :coefname => ByRow(string) => :coefname) + filter!(:coefname => in(_coefnames(x; show_intercept)), df) group = :coefname densvar = :β gdf = groupby(df, group) - y = length(gdf):-1:1 - dens = combine(gdf, densvar => kde => :kde) for (offset, row) in enumerate(reverse(eachrow(dens))) # multiply by 0.95 so that the ridges don't overlap dd = 0.95 * row.kde.density ./ maximum(row.kde.density) - lower = Observable(Point2f.(row.kde.x, offset)) - upper = Observable(Point2f.(row.kde.x, dd .+ offset)) - band!(plot, lower, upper; color=(:black, 0.3)) - lines!(plot, upper; color=(:black, 1.0)) + lower = Point2f.(row.kde.x, offset) + upper = Point2f.(row.kde.x, dd .+ offset) + band!(ax, lower, upper; band_attributes...) + lines!(ax, upper; attributes...) + end + + # check conf_level so that we don't double print + # if coefplot took care of it for us + if ismissing(conf_level) + cn = _coefnames(x; show_intercept) + nticks = _npreds(x; show_intercept) + ax.yticks = (nticks:-1:1, cn) + ylims!(ax, 0, nticks + 1) + vline_at_zero && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) end - return plot + reset_limits!(ax) + + return ax end # """ diff --git a/src/ridge2d.jl b/src/ridge2d.jl index cc15b98..14672c1 100644 --- a/src/ridge2d.jl +++ b/src/ridge2d.jl @@ -1,3 +1,19 @@ +""" + ridge2d(bs::MixedModelBootstrap; ptype=:β, kwargs...) + ridge2d!(f::$(Indexable), bs::MixedModelBootstrap; ptype=:β) + +Plot pairwise bivariate scatter plots with overlain densities for a bootstrap sample. + +`ptype` specifies the set of parameters to examine, e.g. `:β`, `:σ`, `:ρ`. + +The mutating methods return the original object. +""" +function ridge2d(bs::MixedModelBootstrap, args...; kwargs...) + f = Figure(; size=(1000, 1000)) # use an aspect ratio of 1 for the whole figure + + return ridge2d!(f, bs, args...; kwargs...) +end + function _ridge2d_panel!(ax::Axis, i::Int, j::Int, cnames::Vector{Symbol}, tbl) x = Tables.getcolumn(tbl, cnames[j]) y = Tables.getcolumn(tbl, cnames[i]) @@ -10,17 +26,11 @@ function _ridge2d_panel!(ax::Axis, i::Int, j::Int, cnames::Vector{Symbol}, tbl) return plt end -""" - ridge2d!(f::Union{Makie.FigureLike,Makie.GridLayout}, bs::MixedModelBootstrap; - ptype=:β) - -Plot pairwise bivariate scatter plots with overlain densities for a bootstrap sample. +# No Axis variant here because we always need to great a subgrid +# we don't yet expose any further configuration options / attributes - -`ptype` specifies the set of parameters to examine, e.g. `:β`, `:σ`, `:ρ`. -""" -function ridge2d!(f::Union{Makie.FigureLike,Makie.GridLayout}, bs::MixedModelBootstrap; - ptype=:β) +"""$(@doc ridge2d)""" +function ridge2d!(f::Indexable, bs::MixedModelBootstrap; ptype=:β) tbl = bs.tbl cnames = [string(x) for x in propertynames(tbl)[2:end]] filter!(startswith(string(ptype)), cnames) @@ -31,10 +41,3 @@ function ridge2d!(f::Union{Makie.FigureLike,Makie.GridLayout}, bs::MixedModelBoo splomaxes!(f, cnames, _ridge2d_panel!, Symbol.(cnames), tbl) return f end - -"""$(@doc ridge2d!)""" -function ridge2d(bs::MixedModelBootstrap, args...; kwargs...) - f = Figure(; size=(1000, 1000)) # use an aspect ratio of 1 for the whole figure - - return ridge2d!(f, bs, args...; kwargs...) -end diff --git a/src/utilities.jl b/src/utilities.jl index 8ddd398..f532d12 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,16 +1,16 @@ function _coefnames(x; show_intercept=true) - cn = coefnames(x) + cn = fixefnames(x) return show_intercept ? cn : filter!(!=("(Intercept)"), cn) end function _coefnames(x::MixedModelBootstrap; show_intercept=true) - cn = collect(string.(propertynames(first(x.fits).β))) + cn = [string(k) for (k, v) in pairs(first(x.fits).β) if !isequal(v, -0.0)] return show_intercept ? cn : filter!(!=("(Intercept)"), cn) end """ - confint_table(x::StatsBase.StatisticalModel, level=0.95) - confint_table(x::MixedModelBootstrap, level=0.95) + confint_table(x::MixedModel, level=0.95; show_intercept=true) + confint_table(x::MixedModelBootstrap, level=0.95; show_intercept=true) Return a DataFrame of coefficient names, point estimates and confidence intervals. @@ -26,26 +26,33 @@ The returned table has the following columns: - `lower`: the lower edge of the confidence interval - `upper`: the upper edge of the confidence interval +!!! note + This function is internal and may be removed in a future release + without being considered a breaking change. """ -function confint_table(x::StatsBase.StatisticalModel, level=0.95) +function confint_table(x::MixedModel, level=0.95; show_intercept=true) # taking from the lower tail semultiple = zquantile((1 - level) / 2) se = stderror(x) + est = coef(x) - return DataFrame(; - coefname=coefnames(x), - estimate=coef(x), - lower=coef(x) + semultiple * se, - upper=coef(x) - semultiple * se) + df = DataFrame(; + coefname=coefnames(x), + estimate=coef(x), + # signs are 'swapped' b/c semultiple comes from the lower tail + lower=est + semultiple * se, + upper=est - semultiple * se) + return filter!(:coefname => in(_coefnames(x; show_intercept)), df) end -function confint_table(x::MixedModelBootstrap, level=0.95) +function confint_table(x::MixedModelBootstrap, level=0.95; show_intercept=true) df = transform!(select!(DataFrame(x.β), Not(:iter)), :coefname => ByRow(string) => :coefname) ci(x) = shortestcovint(x, level) - return combine(groupby(df, :coefname), - :β => mean => :estimate, - :β => NamedTuple{(:lower, :upper)} ∘ ci => [:lower, :upper]) + df = combine(groupby(df, :coefname), + :β => mean => :estimate, + :β => NamedTuple{(:lower, :upper)} ∘ ci => [:lower, :upper]) + return filter!(:coefname => in(_coefnames(x; show_intercept)), df) end _npreds(x; show_intercept=true) = length(_coefnames(x; show_intercept)) diff --git a/test/Project.toml b/test/Project.toml index cfdfe2a..0400cfa 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,9 +6,11 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04" [compat] Aqua = "0.5, 0.6" +Suppressor = "0.2" TestSetExtensions = "3" diff --git a/test/clevelandaxes.jl b/test/clevelandaxes.jl index 56a1885..0fd1079 100644 --- a/test/clevelandaxes.jl +++ b/test/clevelandaxes.jl @@ -17,6 +17,7 @@ catch err end end for i in 1:4, j in 1:4 + local x, y x = randn(MersenneTwister(i), n) y = randn(MersenneTwister(j), n) scatter!(f[i, j], x, y) diff --git a/test/coefplot.jl b/test/coefplot.jl index 273f1e2..6f5aac6 100644 --- a/test/coefplot.jl +++ b/test/coefplot.jl @@ -3,3 +3,9 @@ f = coefplot(m1) f = coefplot(b1) @test save(joinpath(OUTDIR, "coef_sleepstudy_boot.png"), f) + +f = coefplot(mr) +@test save(joinpath(OUTDIR, "coef_rankdeficient.png"), f) + +f = coefplot(br) +@test save(joinpath(OUTDIR, "coef_rankdeficient_boot.png"), f) diff --git a/test/ridgeplot.jl b/test/ridgeplot.jl index 5cc611a..edc5feb 100644 --- a/test/ridgeplot.jl +++ b/test/ridgeplot.jl @@ -1,2 +1,5 @@ f = ridgeplot(b1) @test save(joinpath(OUTDIR, "ridge_sleepstudy.png"), f) + +f = ridgeplot(br) +@test save(joinpath(OUTDIR, "ridge_rankdeficient.png"), f) diff --git a/test/runtests.jl b/test/runtests.jl index 99c45c6..6831c67 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ using MixedModels using MixedModelsMakie using Random # we don't depend on exact PRNG vals, so no need for StableRNGs using Statistics +using Suppressor using Test using TestSetExtensions @@ -41,6 +42,13 @@ g1 = fit(MixedModel, MixedModels.dataset(:verbagg), Bernoulli(); progress) +rng = MersenneTwister(0) +x = rand(rng, 100) +data = (x=x, x2=1.5 .* x, y=rand(rng, [0, 1], 100), z=repeat('A':'T', 5)) +mr = @suppress fit(MixedModel, @formula(y ~ x + x2 + (1 | z)), data; progress) +br = parametricbootstrap(MersenneTwister(42), 500, mr; progress, + optsum_overrides=(; ftol_rel=1e-6)) + @testset ExtendedTestSet "MixedModelsMakie.jl" begin @testset "Aqua" begin # we can't check for unbound type parameters