diff --git a/src/coefplot.jl b/src/coefplot.jl index 60a02be..8fad660 100644 --- a/src/coefplot.jl +++ b/src/coefplot.jl @@ -1,6 +1,6 @@ """ coefplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure - coefplot!(fig::$(Indexable), x::Union{MixedModel,MixedModelBootstrap}; + 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...) @@ -13,13 +13,13 @@ Inestimable coefficients (coefficients removed by pivoting in the rank deficient The mutating methods return the original object. -!!! note - Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) +!!! note + Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded. """ -function coefplot(x::Union{MixedModel,MixedModelBootstrap}; show_intercept=true, kwargs...) +function coefplot(x::Union{MixedModel,MixedModelBootstrap}; show_intercept=true, ptype=nothing, kwargs...) # need to guarantee a min height of 150 - fig = Figure(; size=(640, max(150, 75 * _npreds(x; show_intercept)))) + fig = Figure(; size=(640, max(150, 75 * _npreds(x, ptype; show_intercept)))) coefplot!(fig, x; kwargs...) return fig end @@ -36,8 +36,9 @@ function coefplot!(ax::Axis, x::Union{MixedModel,MixedModelBootstrap}; conf_level=0.95, vline_at_zero=true, show_intercept=true, + ptype=nothing, attributes...) - ci = confint_table(x, conf_level; show_intercept) + ci = confint_table(x, conf_level; show_intercept, ptype) y = nrow(ci):-1:1 xvals = ci.estimate xlabel = @sprintf "Estimate and %g%% confidence interval" conf_level * 100 @@ -51,8 +52,8 @@ function coefplot!(ax::Axis, x::Union{MixedModel,MixedModelBootstrap}; vline_at_zero && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) reset_limits!(ax) - cn = _coefnames(x; show_intercept) - nticks = _npreds(x; show_intercept) + cn = _coefnames(x, ptype; show_intercept) + nticks = _npreds(x, ptype; show_intercept) ax.yticks = (nticks:-1:1, cn) ylims!(ax, 0, nticks + 1) return ax diff --git a/src/ridge.jl b/src/ridge.jl index d2ad163..e3f4803 100644 --- a/src/ridge.jl +++ b/src/ridge.jl @@ -1,6 +1,6 @@ """ ridgeplot(x::Union{MixedModel,MixedModelBootstrap}; kwargs...)::Figure - ridgeplot!(fig::$(Indexable), x::Union{MixedModel,MixedModelBootstrap}; + 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...) @@ -16,13 +16,13 @@ Setting `conf_level=missing` removes the markings for the highest density interv The mutating methods return the original object. -!!! note - Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) +!!! note + Inestimable coefficients (coefficients removed by pivoting in the rank deficient case) are excluded. """ -function ridgeplot(x::MixedModelBootstrap; show_intercept=true, kwargs...) +function ridgeplot(x::MixedModelBootstrap; show_intercept=true, ptype=nothing, kwargs...) # need to guarantee a min height of 200 - fig = Figure(; size=(640, max(200, 100 * _npreds(x; show_intercept)))) + fig = Figure(; size=(640, max(200, 100 * _npreds(x, ptype; show_intercept)))) return ridgeplot!(fig, x; show_intercept, kwargs...) end @@ -36,7 +36,7 @@ end """ _color(s::Symbol) _color(p::Pair) - + Extract the color part out of either a color name or a `(color, alpha)` pair. """ _color(s) = s @@ -47,6 +47,7 @@ function ridgeplot!(ax::Axis, x::MixedModelBootstrap; conf_level=0.95, vline_at_zero=true, show_intercept=true, + ptype=:β, attributes...) xlabel = if !ismissing(conf_level) @sprintf "Normalized bootstrap density and %g%% confidence interval" (conf_level * @@ -56,19 +57,22 @@ function ridgeplot!(ax::Axis, x::MixedModelBootstrap; end if !ismissing(conf_level) - coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, color=:black, - attributes...) + coefplot!(ax, x; conf_level, vline_at_zero, show_intercept, ptype, + color=:black, attributes...) end 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) + df = DataFrame(getproperty(boot, ptype)) + rename!(c -> replace(c, "column" => "coefname"), df) + transform!(df, :coefname => ByRow(string) => :coefname) filter!(:coefname => in(_coefnames(x; show_intercept)), df) group = :coefname - densvar = :β + # drop trailing s + densvar = replace(string(ptype), "s" => "") + # TODO: coefplot gdf = groupby(df, group) dens = combine(gdf, densvar => kde => :kde) @@ -96,75 +100,3 @@ function ridgeplot!(ax::Axis, x::MixedModelBootstrap; return ax end - -# """ -# ridgeplot!(ax::Axis, df::AbstractDataFrame, densvar::Symbol, group::Symbol; normalize=false) -# ridgeplot!(f::Union{Makie.FigureLike,Makie.GridLayout}, args...; pos=(1,1) kwargs...) -# ridgeplot(args...; kwargs...) - -# Create a "ridge plot". - -# A ridge plot is stacked plot of densities for a given variable (`densvar`) grouped by a different variable (`group`). Because densities can very widely in scale, it is sometimes useful to `normalize` the densities so that each density has a maximum of 1. - -# The non-mutating method creates a Figure before calling the method for Figure. -# The method for Figure places the ridge plot in the grid position specified by `pos`, default is (1,1). -# """ -# function ridgeplot!(ax::Axis, df::AbstractDataFrame, densvar::Symbol, group::Symbol; sort_by_group=false, vline_at_zero=true, normalize=false) -# # `normalize` makes it so that the max density is always 1 -# # `normalize` works on the density not the area/mass -# gdf = groupby(df, group) -# dens = combine(gdf, densvar => kde => :kde) -# sort_by_group && sort!(dens, group) -# spacing = normalize ? 1.0 : 0.9 * maximum(dens[!, :kde]) do val -# return maximum(val.density) -# end - -# nticks = length(gdf) - -# for (idx, row) in enumerate(eachrow(dens)) -# dd = normalize ? row.kde.density ./ maximum(row.kde.density) : row.kde.density - -# offset = idx * spacing - -# lower = Node(Point2f.(row.kde.x, offset)) -# upper = Node(Point2f.(row.kde.x, dd .+ offset)) -# band!(ax, lower, upper; color=(:black, 0.3)) -# lines!(ax, upper; color=(:black, 1.0)) -# end - -# vline_at_zero && vlines!(ax, 0; color=(:black, 0.75), linestyle=:dash) - -# ax.yticks[] = (1:spacing:(nticks*spacing), string.(dens[!, group])) -# ylims!(ax, 0, (nticks + 2) * spacing) -# ax.xlabel[] = string(densvar) -# ax.ylabel[] = string(group) - -# ax -# end - -# function ridgeplot!(f::Union{Makie.FigureLike,Makie.GridLayout}, args...; pos=(1,1), kwargs...) -# ridgeplot!(Axis(f[pos...]), args...; kwargs...) -# return f -# end - -# """ -# ridgeplot(args...; kwargs...) - -# See [ridgeplot!](@ref). -# """ -# ridgeplot(args...; kwargs...) = ridgeplot!(Figure(), args...; kwargs...) - -# """ -# ridgeplot!(Union{Figure, Axis}, bstrp::MixedModelBootstrap, args...; show_intercept=true, kwargs...) -# ridgeplot(bstrp::MixedModelBootstrap, args...; show_intercept=true, kwargs...) - -# Convenience methods that call `DataFrame(bstrp.β)` and pass that onto the primary `ridgeplot[!]` methods. -# By default, the intercept is shown, but this can be disabled. -# """ -# function ridgeplot!(axis::Axis, bstrp::MixedModelBootstrap, args...; -# show_intercept=true, normalize=true, kwargs...) - -# df = DataFrame(bstrp.β) -# show_intercept || filter!(:coefname => !=(Symbol("(Intercept)")), df) -# ridgeplot!(axis, df, :β, :coefname, args...; sort_by_group=false, normalize, kwargs...) -# end diff --git a/src/utilities.jl b/src/utilities.jl index f532d12..59341b8 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,10 +1,12 @@ -function _coefnames(x; show_intercept=true) +function _coefnames(x::MixedModel, ptype::Nothing=nothing; show_intercept=true) cn = fixefnames(x) return show_intercept ? cn : filter!(!=("(Intercept)"), cn) end -function _coefnames(x::MixedModelBootstrap; show_intercept=true) - cn = [string(k) for (k, v) in pairs(first(x.fits).β) if !isequal(v, -0.0)] +function _coefnames(x::MixedModelBootstrap, ptype; show_intercept=true) + ptype = something(ptype, :β) + nt = getproperty(first(x.fits), ptype) + cn = [string(k) for (k, v) in pairs(nt) if !isequal(v, -0.0)] return show_intercept ? cn : filter!(!=("(Intercept)"), cn) end @@ -30,7 +32,8 @@ The returned table has the following columns: This function is internal and may be removed in a future release without being considered a breaking change. """ -function confint_table(x::MixedModel, level=0.95; show_intercept=true) +function confint_table(x::MixedModel, level=0.95; ptype=nothing, show_intercept=true) + isnothing(ptype) || throw(ArgumentError("`ptype` not supported for MixedModel")) # taking from the lower tail semultiple = zquantile((1 - level) / 2) se = stderror(x) @@ -45,17 +48,23 @@ function confint_table(x::MixedModel, level=0.95; show_intercept=true) return filter!(:coefname => in(_coefnames(x; show_intercept)), df) end -function confint_table(x::MixedModelBootstrap, level=0.95; show_intercept=true) - df = transform!(select!(DataFrame(x.β), Not(:iter)), - :coefname => ByRow(string) => :coefname) +function confint_table(x::MixedModelBootstrap, level=0.95; ptype=:β, show_intercept=true) + ptype = something(ptype, :β) + ptype in (:β, :σs) || throw(ArgumentError("ptype $(ptype) not supported")) + df = DataFrame(getproperty(x, ptype)) + rename!(c -> replace(c, "column" => "coefname"), df) + transform!(select!(df, Not(:iter)), + :coefname => ByRow(string) => :coefname) ci(x) = shortestcovint(x, level) + # drop trailing s + var = replace(string(ptype), "s" => "") df = combine(groupby(df, :coefname), - :β => mean => :estimate, - :β => NamedTuple{(:lower, :upper)} ∘ ci => [:lower, :upper]) - return filter!(:coefname => in(_coefnames(x; show_intercept)), df) + var => mean => :estimate, + var => NamedTuple{(:lower, :upper)} ∘ ci => [:lower, :upper]) + return filter!(:coefname => in(_coefnames(x, ptype; show_intercept)), df) end -_npreds(x; show_intercept=true) = length(_coefnames(x; show_intercept)) +_npreds(args...; kwargs...) = length(_coefnames(args...; kwargs...)) """ ppoints(n::Integer)