From e701d8398c6c2969d913a79f46c235efb85e7941 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Sun, 8 Dec 2024 22:42:03 +0100 Subject: [PATCH 01/13] Plots stub --- ext/PlotsExt.jl | 67 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 ext/PlotsExt.jl diff --git a/ext/PlotsExt.jl b/ext/PlotsExt.jl new file mode 100644 index 00000000..aba15d40 --- /dev/null +++ b/ext/PlotsExt.jl @@ -0,0 +1,67 @@ +module PlotsExt + using Plots + using RecipesBase + + @userplot struct LMPlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(lmp::LMPlot) + + end + + @userplot struct ResidualPlot{T<:Tuple} + args::T + end + + @recipe function f(rp::ResidualPlot) + xlabel --> "Fitted values" + ylabel --> "Residuals" + title --> "Residuals vs Fitted" + end + + @userplot struct ScaleLocationPlot{T<:Tuple} + args::T + end + + @recipe function f(slp::ScaleLocationPlot) + xlabel --> "Fitted values" + ylabel --> L"\sqrt{|\text{standardized residuals}|}" + title --> "Scale-Location" + end + + @userplot struct QQPlot{T<:Tuple} + args::T + end + + @recipe function f(qqp::QQPlot; qqline=true) + xlabel --> "Theoretical Quantiles" + ylabel --> "Standardized residuals" + title --> "Q-Q Residuals" + + r = residuals(qqp.args[1]) + end + + + @userplot struct ResidualsLeveragePlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(rlp::ResidualsLeveragePlot) + xlabel --> "Leverage" + ylabel --> "Standardized residuals" + title --> "Residuals vs Leverage" + + r = residuals(rlp.args[1]) + end + + @userplot struct CooksLeveragePlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(clp::CooksLeveragePlot) + xlabel --> L"\text{Leverage } h_{ii}" + ylabel --> "Cook's distance" + title --> L"\text{Cook's distance vs Leverage }*h_{ii}/(1-h_{ii})" + end +end From 45c5cd914fd0afb8716ea4e199f706738575cce6 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 24 Dec 2024 23:52:34 +0100 Subject: [PATCH 02/13] Add Basic Plots Recipes --- Project.toml | 10 ++ ext/MakieExt.jl | 3 + ext/StatsPlotsExt.jl | 213 +++++++++++++++++++++++++++++++++++++++++++ src/GLM.jl | 15 +++ src/lm.jl | 13 ++- 5 files changed, 251 insertions(+), 3 deletions(-) create mode 100644 ext/MakieExt.jl create mode 100644 ext/StatsPlotsExt.jl diff --git a/Project.toml b/Project.toml index ddd8ca51..d3d93c31 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "2.0.0-DEV" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -16,6 +17,14 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +[weakdeps] +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" + +[extensions] +StatsPlotsExt = "StatsPlots" +MakieExt = "Makie" + [compat] CSV = "0.7, 0.8, 0.9, 0.10" CategoricalArrays = "0.8, 0.9, 0.10" @@ -29,6 +38,7 @@ StatsAPI = "1.4" StatsBase = "0.33.5, 0.34" StatsFuns = "0.6, 0.7, 0.8, 0.9, 1.0" StatsModels = "0.7.3" +StatsPlots = "0.15" Tables = "1" julia = "1.6" diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl new file mode 100644 index 00000000..ae14bddb --- /dev/null +++ b/ext/MakieExt.jl @@ -0,0 +1,3 @@ +module MakieExt + +end diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl new file mode 100644 index 00000000..71f664d8 --- /dev/null +++ b/ext/StatsPlotsExt.jl @@ -0,0 +1,213 @@ +module StatsPlotsExt + using GLM + using Statistics + using StatsPlots + using RecipesBase + using Distributions + using GLM: leverage + import GLM: cooksleverageplot, cooksleverageplot! + import GLM: scalelocationplot, scalelocationplot! + import GLM: residualplot, residualplot! + import GLM: residualsleverageplot, residualsleverageplot! + import StatsPlots: QQPlot, QQNorm + import StatsPlots: qqplot, qqplot!, qqnorm, qqnorm! + + + function standardized_residuals(obj::LinearModel) + r = residuals(obj) + h = leverage(obj) + return r ./(std(r) .* sqrt.(1 .- h)) + end + + @recipe function f(l::LinearModel) + + end + + @userplot struct ResidualPlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(rp::ResidualPlot; ) + xlabel --> "Fitted values" + ylabel --> "Residuals" + title --> "Residuals vs Fitted" + label --> "" + + r = residuals(rp.args[1]) + y = predict(rp.args[1]) + + @series begin + seriestype := :scatter + y, r + end + + @series begin + seriestype := :hline + linecolor := :black + linestyle := :dash + linewidth := 0.5 + label := "" + [0.0] + end + nothing + end + + @userplot struct ScaleLocationPlot{T<:Tuple} + args::T + end + + @recipe function f(slp::ScaleLocationPlot) + xlabel --> "Fitted values" + ylabel --> "√|standardized residuals|" + title --> "Scale-Location" + label --> "" + + r = standardized_residuals(slp.args[1]) + y = predict(slp.args[1]) + @series begin + seriestype := :scatter + y, (sqrt ∘ abs).(r) + end + nothing + end + + #= + @userplot struct QuantileQuantilePlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(qqp::QuantileQuantilePlot) + xlabel --> "Theoretical Quantiles" + ylabel --> "Standardized residuals" + title --> "Q-Q Residuals" + label --> "" + + r = residuals(qqp.args[1]) + @series begin + seriestype := :qqnorm + linestyle := :dash + linecolor := :black + linewidth := 0.5 + r + end + end + =# + @recipe function f(::QQPlot, obj::LinearModel) + xlabel --> "Theoretical Quantiles" + ylabel --> "Standardized Residuals" + title --> "Q-Q Residuals" + label --> "" + + linestyle --> :dash + linecolor --> :gray + linewidth --> 0.5 + + QQPlot(Normal, standardized_residuals(obj)) + end + + # This feels hacky but it works + qqplot(l::LinearModel, D=Normal; + xlabel = "Theoretical Quantiles", + ylabel = "Standardized Residuals", + title = "Q-Q Residuals", + label = "", + linestyle = :dash, + linecolor = :gray, + linewidth = 0.5, + kw... + ) = qqplot(D, standardized_residuals(l); lsty=linestyle, lcol=linecolor, lw=linewidth, xlabel=xlabel,ylabel=ylabel,title=title,label=label, kw...) + + qqplot!(p, l::LinearModel, D=Normal; kw...) = qqplot!(p, D, standardized_residuals(l); kw...) + + qqnorm(l::LinearModel; kw...) = qqnorm(standardized_residuals(l); kw...) + qqnorm!(p, l::LinearModel; kw...) = qqnorm!(p, standardized_residuals(l), kw...) + + + @userplot struct ResidualsLeveragePlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(rlp::ResidualsLeveragePlot; cook_levels = [0.5, 1.0]) + xlabel --> "Leverage" + ylabel --> "Standardized residuals" + title --> "Residuals vs Leverage" + label --> "" + + r = standardized_residuals(rlp.args[1]) + h = leverage(rlp.args[1]) + k = dof(rlp.args[1]) - 1 + + ymax = maximum(abs.(r)) + ylims --> (-1.1*ymax, 1.1*ymax) + + @series begin + seriestype := :scatter + h, r + end + + @series begin + seriestype := :hline + linecolor := :black + linestyle := :dash + linewidth := 0.5 + label := "" + [0.0] + end + @series begin + seriestype := :vline + linecolor := :black + linestyle := :dash + linewidth := 0.5 + [0.0] + end + if !isempty(cook_levels) + cookfun = (h,D,k) -> sqrt(D * k * (1-h) / h) + xmin, xmax = extrema(h) + xs = LinRange(xmin, 1.1*xmax, 50) + for D in cook_levels + @series begin + seriestype := :path + linecolor := :gray + linestyle := :dash + annotations := [(xs[end], cookfun(xs[end],D,k), ("$D", 8, :gray))] + xs, cookfun.(xs,D,k) + end + + @series begin + seriestype := :scatter + markeralpha := 0.0 + series_annotation := [("$D", 9, :gray)] + [1.02 * xs[end]], [cookfun(xs[end],D,k)] + end + + @series begin + seriestype := :path + linecolor := :gray + linestyle := :dash + xs, -cookfun.(xs,D,k) + end + end + end + + + end + + @userplot struct CooksLeveragePlot{T<:Tuple{LinearModel}} + args::T + end + + @recipe function f(clp::CooksLeveragePlot) + xlabel --> "Leverage" + ylabel --> "Cook's distance" + title --> "Cook's distance vs Leverage" + label --> "" + + + h = leverage(clp.args[1]) + cd = cooksdistance(clp.args[1]) + @series begin + seriestype := :scatter + h, cd + end + end +end diff --git a/src/GLM.jl b/src/GLM.jl index 59c327db..c218c801 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -67,6 +67,21 @@ module GLM predict, # make predictions ftest # compare models with an F test + # Plot functions + export cooksleverageplot, cooksleverageplot! + export scalelocationplot, scalelocationplot! + export residualplot, residualplot! + export residualsleverageplot, residualsleverageplot! + export quantilequantileqplot, quantilequantileplot! + function cooksleverageplot end + function cooksleverageplot! end + function scalelocationplot end + function scalelocationplot! end + function residualplot end + function residualplot! end + function residualsleverageplot end + function residualsleverageplot! end + const FP = AbstractFloat const FPVector{T<:FP} = AbstractArray{T,1} diff --git a/src/lm.jl b/src/lm.jl index f26fb8e5..a4f085e9 100644 --- a/src/lm.jl +++ b/src/lm.jl @@ -353,14 +353,21 @@ function StatsBase.cooksdistance(obj::LinearModel) k = dof(obj)-1 d_res = dof_residual(obj) X = modelmatrix(obj) - XtX = crossmodelmatrix(obj) k == size(X,2) || throw(ArgumentError("Models with collinear terms are not currently supported.")) + hii = leverage(obj) + + D = @. u^2 * (hii / (1 - hii)^2) / (k*mse) + return D +end + +function leverage(obj::LinearModel) + X = modelmatrix(obj) + XtX = crossmodelmatrix(obj) wts = obj.rr.wts if isempty(wts) hii = diag(X * inv(XtX) * X') else throw(ArgumentError("Weighted models are not currently supported.")) end - D = @. u^2 * (hii / (1 - hii)^2) / (k*mse) - return D + return hii end From c2c9ad0219f9358e1e4162b3f67896ce60fa3066 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Fri, 27 Dec 2024 00:00:27 +0100 Subject: [PATCH 03/13] Add Makie recipes --- ext/MakieExt.jl | 174 +++++++++++++++++++++++++++++++++++++++++++ ext/PlotsExt.jl | 67 ----------------- ext/StatsPlotsExt.jl | 26 ++++--- src/GLM.jl | 10 +-- src/plots.jl | 26 +++++++ 5 files changed, 216 insertions(+), 87 deletions(-) delete mode 100644 ext/PlotsExt.jl create mode 100644 src/plots.jl diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl index ae14bddb..46116fea 100644 --- a/ext/MakieExt.jl +++ b/ext/MakieExt.jl @@ -1,3 +1,177 @@ module MakieExt + using Makie + using GLM + using Distributions + using GLM: leverage, standardized_residuals + import GLM: cooksleverageplot, cooksleverageplot! + import GLM: scalelocationplot, scalelocationplot! + import GLM: residualplot, residualplot! + import GLM: residualsleverageplot, residualsleverageplot! + import GLM: quantilequantileplot, quantilequantileplot! + Makie.@recipe(ResidualPlot, obj) do scene + Theme( + axislines = true, + axislinestyle = :dot, + axislinecolor = :gray, + axislinewidth = 1, + axislabels = true, + axistitle = true + ) + #Attributes() + end + + function Makie.plot!(rp::ResidualPlot{<:Tuple{LinearModel}}) + r = residuals(rp.obj[]) + y = predict(rp.obj[]) + + ax = current_axis() + + if rp.axistitle[] + ax.title = "Residuals vs Fitted" + end + if rp.axislabels[] + ax.xlabel = "Fitted Values" + ax.ylabel = "Residuals" + end + + scatter!(rp, y, r) + if rp.axislines[] + hlines!(rp, [0], color = rp.axislinecolor[], linestyle = rp.axislinestyle[], linewidth = rp.axislinewidth[]) + end + rp + end + + Makie.@recipe(ScaleLocationPlot, obj) do scene + Theme( + axislabels = true, + axistitle = true + ) + end + + function Makie.plot!(slp::ScaleLocationPlot{<:Tuple{LinearModel}}) + ax = current_axis() + if slp.axistitle[] + ax.title = "Scale-Location" + end + if slp.axislabels[] + ax.xlabel = "Fitted Values" + ax.ylabel = "√|standardized residuals|" + end + + r = standardized_residuals(slp.obj[]) + y = predict(slp.obj[]) + + scatter!(slp, y, (sqrt ∘ abs).(r)) + + return slp + end + + Makie.@recipe(ResidualsLeveragePlot, obj) do scene + Theme( + axislines = true, + axislinestyle = :dot, + axislinecolor = :gray, + axislinewidth = 1, + axislabels = true, + axistitle = true, + cookslevels = [0.5,1.0], + cookslinecolor = :gray, + cookslinestyle = :dash, + cookslinewidth = 1 + ) + end + + function Makie.plot!(rlp::ResidualsLeveragePlot{<:Tuple{LinearModel}}) + ax = current_axis() + + if rlp.axistitle[] + ax.title = "Residuals vs Leverage" + end + if rlp.axislabels[] + ax.xlabel = "Leverage" + ax.ylabel = "Standardized Residuals" + end + + r = standardized_residuals(rlp.obj[]) + h = leverage(rlp.obj[]) + k = dof(rlp.obj[]) - 1 + + ymax = maximum(abs.(r)) + ylims!(ax, -1.1*ymax, 1.1*ymax) + + if rlp.axislines[] + hlines!(rlp, [0], + linestyle = rlp.axislinestyle[], + color = rlp.axislinecolor[], + linewidth = rlp.axislinewidth[] + ) + vlines!(rlp, [0], + linestyle = rlp.axislinestyle[], + color = rlp.axislinecolor[], + linewidth = rlp.axislinewidth[] + ) + end + + if !isempty(rlp.cookslevels[]) + cooksfun = (h,D,k) -> sqrt(D * k * (1-h) / h) + xmin, xmax = extrema(h) + xs = LinRange(xmin, 1.1*xmax, 50) + for D in rlp.cookslevels[] + lines!(rlp, xs, cooksfun.(xs,D,k), + linestyle = rlp.cookslinestyle[], + color = rlp.cookslinecolor[], + linewidth = rlp.cookslinewidth[] + ) + lines!(rlp, xs, -cooksfun.(xs,D,k), + linestyle = rlp.cookslinestyle[], + color = rlp.cookslinecolor[], + linewidth = rlp.cookslinewidth[] + ) + end + end + + scatter!(rlp, h, r) + return rlp + end + + Makie.@recipe(CooksLeveragePlot, obj) do scene + Theme( + axistitle = true, + axislabels = true, + ) + end + + function Makie.plot!(clp::CooksLeveragePlot{<:Tuple{LinearModel}}) + ax = current_axis() + if clp.axistitle[] + ax.title = "Cook's distance vs Leverage" + end + if clp.axislabels[] + ax.xlabel = "Leverage" + ax.ylabel = "Cook's Distance" + end + + h = leverage(clp.obj[]) + cd = cooksdistance(clp.obj[]) + + scatter!(clp, h,cd) + return clp + end + + function Makie.plot!(qqp::QQPlot{<:Tuple{LinearModel}}) + ax = current_axis() + + ax.ylabel = "Standardized Residuals" + ax.xlabel = "Theoretical Quantiles" + ax.title = "Q-Q Residuals" + + r = standardized_residuals(qqp[1][]) + qqplot!(qqp, Normal, r, + qqline = :identity, + linestyle = :dash, + linewidth = 1 + ) + return qqp + end end diff --git a/ext/PlotsExt.jl b/ext/PlotsExt.jl deleted file mode 100644 index aba15d40..00000000 --- a/ext/PlotsExt.jl +++ /dev/null @@ -1,67 +0,0 @@ -module PlotsExt - using Plots - using RecipesBase - - @userplot struct LMPlot{T<:Tuple{LinearModel}} - args::T - end - - @recipe function f(lmp::LMPlot) - - end - - @userplot struct ResidualPlot{T<:Tuple} - args::T - end - - @recipe function f(rp::ResidualPlot) - xlabel --> "Fitted values" - ylabel --> "Residuals" - title --> "Residuals vs Fitted" - end - - @userplot struct ScaleLocationPlot{T<:Tuple} - args::T - end - - @recipe function f(slp::ScaleLocationPlot) - xlabel --> "Fitted values" - ylabel --> L"\sqrt{|\text{standardized residuals}|}" - title --> "Scale-Location" - end - - @userplot struct QQPlot{T<:Tuple} - args::T - end - - @recipe function f(qqp::QQPlot; qqline=true) - xlabel --> "Theoretical Quantiles" - ylabel --> "Standardized residuals" - title --> "Q-Q Residuals" - - r = residuals(qqp.args[1]) - end - - - @userplot struct ResidualsLeveragePlot{T<:Tuple{LinearModel}} - args::T - end - - @recipe function f(rlp::ResidualsLeveragePlot) - xlabel --> "Leverage" - ylabel --> "Standardized residuals" - title --> "Residuals vs Leverage" - - r = residuals(rlp.args[1]) - end - - @userplot struct CooksLeveragePlot{T<:Tuple{LinearModel}} - args::T - end - - @recipe function f(clp::CooksLeveragePlot) - xlabel --> L"\text{Leverage } h_{ii}" - ylabel --> "Cook's distance" - title --> L"\text{Cook's distance vs Leverage }*h_{ii}/(1-h_{ii})" - end -end diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 71f664d8..791f946c 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -4,20 +4,18 @@ module StatsPlotsExt using StatsPlots using RecipesBase using Distributions - using GLM: leverage + using GLM: leverage, standardized_residuals + using RecipesBase: recipetype import GLM: cooksleverageplot, cooksleverageplot! import GLM: scalelocationplot, scalelocationplot! import GLM: residualplot, residualplot! import GLM: residualsleverageplot, residualsleverageplot! + import GLM: quantilequantileplot, quantilequantileplot! import StatsPlots: QQPlot, QQNorm import StatsPlots: qqplot, qqplot!, qqnorm, qqnorm! - function standardized_residuals(obj::LinearModel) - r = residuals(obj) - h = leverage(obj) - return r ./(std(r) .* sqrt.(1 .- h)) - end + @recipe function f(l::LinearModel) @@ -52,7 +50,7 @@ module StatsPlotsExt nothing end - @userplot struct ScaleLocationPlot{T<:Tuple} + @userplot struct ScaleLocationPlot{T<:Tuple{LinearModel}} args::T end @@ -71,11 +69,11 @@ module StatsPlotsExt nothing end - #= + @userplot struct QuantileQuantilePlot{T<:Tuple{LinearModel}} args::T end - + #= @recipe function f(qqp::QuantileQuantilePlot) xlabel --> "Theoretical Quantiles" ylabel --> "Standardized residuals" @@ -92,7 +90,8 @@ module StatsPlotsExt end end =# - @recipe function f(::QQPlot, obj::LinearModel) + + @recipe function f(qqp::QuantileQuantilePlot) xlabel --> "Theoretical Quantiles" ylabel --> "Standardized Residuals" title --> "Q-Q Residuals" @@ -102,9 +101,12 @@ module StatsPlotsExt linecolor --> :gray linewidth --> 0.5 - QQPlot(Normal, standardized_residuals(obj)) + QQPlot((Normal, standardized_residuals(qqp.args[1]))) end + StatsPlots.recipetype(::Val{:qqplot}, obj::LinearModel, args...) = QuantileQuantilePlot((obj, args...)) + + #= # This feels hacky but it works qqplot(l::LinearModel, D=Normal; xlabel = "Theoretical Quantiles", @@ -121,7 +123,7 @@ module StatsPlotsExt qqnorm(l::LinearModel; kw...) = qqnorm(standardized_residuals(l); kw...) qqnorm!(p, l::LinearModel; kw...) = qqnorm!(p, standardized_residuals(l), kw...) - + =# @userplot struct ResidualsLeveragePlot{T<:Tuple{LinearModel}} args::T diff --git a/src/GLM.jl b/src/GLM.jl index c218c801..7a8b786d 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -73,14 +73,7 @@ module GLM export residualplot, residualplot! export residualsleverageplot, residualsleverageplot! export quantilequantileqplot, quantilequantileplot! - function cooksleverageplot end - function cooksleverageplot! end - function scalelocationplot end - function scalelocationplot! end - function residualplot end - function residualplot! end - function residualsleverageplot end - function residualsleverageplot! end + const FP = AbstractFloat const FPVector{T<:FP} = AbstractArray{T,1} @@ -147,5 +140,6 @@ module GLM include("ftest.jl") include("negbinfit.jl") include("deprecated.jl") + include("plots.jl") end # module diff --git a/src/plots.jl b/src/plots.jl new file mode 100644 index 00000000..bb3a1c36 --- /dev/null +++ b/src/plots.jl @@ -0,0 +1,26 @@ + +""" + standardized_residuals(obj::LinearModel) + +Compute the standardized residuals of a linear model, defined for the `i`-th observation as +``` +r[i] / (std(r) * sqrt(1 - h[i]), +``` +where `r` are the residuals of the model, ``s`` is the empirical standard deviation of the residuals and ``h[i]`` is the leverage of observation `i`. +""" +function standardized_residuals(obj::LinearModel) + r = residuals(obj) + h = leverage(obj) + return r ./(std(r) .* sqrt.(1 .- h)) +end + +function cooksleverageplot end +function cooksleverageplot! end +function scalelocationplot end +function scalelocationplot! end +function residualplot end +function residualplot! end +function residualsleverageplot end +function residualsleverageplot! end +function quantilequantileplot end +function quantilequantileplot! end From 9450a047913a63aa7d43d3f15ccb2d8bfb18037e Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 31 Dec 2024 18:00:33 +0100 Subject: [PATCH 04/13] Makie Recipes + lmplot --- ext/MakieExt.jl | 126 ++++++++++++++++++++++++------------------------ src/GLM.jl | 2 +- src/plots.jl | 2 + 3 files changed, 67 insertions(+), 63 deletions(-) diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl index 46116fea..0b2cc129 100644 --- a/ext/MakieExt.jl +++ b/ext/MakieExt.jl @@ -7,7 +7,51 @@ module MakieExt import GLM: scalelocationplot, scalelocationplot! import GLM: residualplot, residualplot! import GLM: residualsleverageplot, residualsleverageplot! - import GLM: quantilequantileplot, quantilequantileplot! + import GLM: lmplot + import Makie: qqplot, qqplot! + + + function lmplot(obj::LinearModel; figkw...) + fig = Figure(; figkw...) + ax_1 = Axis(fig[1,1], + title = "Residuals vs Fitted Values", + xlabel = "Fitted Values", + ylabel = "Residuals" + ) + residualplot!(ax_1, obj) + ax_2 = Axis(fig[1,2], + title = "Q-Q Residuals", + xlabel = "Theoretical Quantiles", + ylabel = "Standardized Residuals" + ) + qqplot!(ax_2, obj) + ax_3 = Axis(fig[2,1], + title = "Scale-Location", + xlabel = "Fitted Values", + ylabel = "√|standardized residuals|" + ) + scalelocationplot!(ax_3, obj) + ax_4 = Axis(fig[2,2], + title = "Residuals vs Leverage", + xlabel = "Leverage", + ylabel = "Standardized Residuals" + ) + ymax = maximum(abs.(standardized_residuals(obj))) + xmax = maximum(leverage(obj)) + ylims!(ax_4, -1.2*ymax, 1.2*ymax) + xlims!(ax_4, 0.0, 1.2*xmax) + axislegend(ax_4, + [LineElement(color = :gray, linestyle = :dash, linewidth=1)], + ["Cook's distance"], + position = :lb, + framevisible = false, + labelsize = 10, + labelcolor = :gray, + padding = (0.0f0, 0.0f0, 0.0f0, 0.0f0) + ) + residualsleverageplot!(ax_4, obj) + fig + end Makie.@recipe(ResidualPlot, obj) do scene Theme( @@ -15,26 +59,14 @@ module MakieExt axislinestyle = :dot, axislinecolor = :gray, axislinewidth = 1, - axislabels = true, - axistitle = true ) - #Attributes() + end function Makie.plot!(rp::ResidualPlot{<:Tuple{LinearModel}}) r = residuals(rp.obj[]) y = predict(rp.obj[]) - ax = current_axis() - - if rp.axistitle[] - ax.title = "Residuals vs Fitted" - end - if rp.axislabels[] - ax.xlabel = "Fitted Values" - ax.ylabel = "Residuals" - end - scatter!(rp, y, r) if rp.axislines[] hlines!(rp, [0], color = rp.axislinecolor[], linestyle = rp.axislinestyle[], linewidth = rp.axislinewidth[]) @@ -43,27 +75,14 @@ module MakieExt end Makie.@recipe(ScaleLocationPlot, obj) do scene - Theme( - axislabels = true, - axistitle = true - ) + Theme() end function Makie.plot!(slp::ScaleLocationPlot{<:Tuple{LinearModel}}) - ax = current_axis() - if slp.axistitle[] - ax.title = "Scale-Location" - end - if slp.axislabels[] - ax.xlabel = "Fitted Values" - ax.ylabel = "√|standardized residuals|" - end - r = standardized_residuals(slp.obj[]) y = predict(slp.obj[]) scatter!(slp, y, (sqrt ∘ abs).(r)) - return slp end @@ -73,8 +92,6 @@ module MakieExt axislinestyle = :dot, axislinecolor = :gray, axislinewidth = 1, - axislabels = true, - axistitle = true, cookslevels = [0.5,1.0], cookslinecolor = :gray, cookslinestyle = :dash, @@ -83,22 +100,11 @@ module MakieExt end function Makie.plot!(rlp::ResidualsLeveragePlot{<:Tuple{LinearModel}}) - ax = current_axis() - - if rlp.axistitle[] - ax.title = "Residuals vs Leverage" - end - if rlp.axislabels[] - ax.xlabel = "Leverage" - ax.ylabel = "Standardized Residuals" - end - r = standardized_residuals(rlp.obj[]) h = leverage(rlp.obj[]) k = dof(rlp.obj[]) - 1 - ymax = maximum(abs.(r)) - ylims!(ax, -1.1*ymax, 1.1*ymax) + scatter!(rlp, h, r) if rlp.axislines[] hlines!(rlp, [0], @@ -123,35 +129,36 @@ module MakieExt color = rlp.cookslinecolor[], linewidth = rlp.cookslinewidth[] ) + text!(rlp, xs[end], cooksfun(xs[end],D,k), + text = "$D", + color = rlp.cookslinecolor[], + #align = (:left, :bottom), + offset = (1,-4), + fontsize = 10 + ) lines!(rlp, xs, -cooksfun.(xs,D,k), linestyle = rlp.cookslinestyle[], color = rlp.cookslinecolor[], linewidth = rlp.cookslinewidth[] ) + text!(rlp, xs[end], -cooksfun(xs[end],D,k), + text = "$D", + color = rlp.cookslinecolor[], + #align = (:left, :bottom), + offset = (1,-4), + fontsize = 10 + ) end end - scatter!(rlp, h, r) return rlp end Makie.@recipe(CooksLeveragePlot, obj) do scene - Theme( - axistitle = true, - axislabels = true, - ) + Theme() end function Makie.plot!(clp::CooksLeveragePlot{<:Tuple{LinearModel}}) - ax = current_axis() - if clp.axistitle[] - ax.title = "Cook's distance vs Leverage" - end - if clp.axislabels[] - ax.xlabel = "Leverage" - ax.ylabel = "Cook's Distance" - end - h = leverage(clp.obj[]) cd = cooksdistance(clp.obj[]) @@ -160,12 +167,6 @@ module MakieExt end function Makie.plot!(qqp::QQPlot{<:Tuple{LinearModel}}) - ax = current_axis() - - ax.ylabel = "Standardized Residuals" - ax.xlabel = "Theoretical Quantiles" - ax.title = "Q-Q Residuals" - r = standardized_residuals(qqp[1][]) qqplot!(qqp, Normal, r, qqline = :identity, @@ -174,4 +175,5 @@ module MakieExt ) return qqp end + end diff --git a/src/GLM.jl b/src/GLM.jl index 7a8b786d..72502dae 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -72,7 +72,7 @@ module GLM export scalelocationplot, scalelocationplot! export residualplot, residualplot! export residualsleverageplot, residualsleverageplot! - export quantilequantileqplot, quantilequantileplot! + export lmplot const FP = AbstractFloat diff --git a/src/plots.jl b/src/plots.jl index bb3a1c36..29eece21 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -14,6 +14,8 @@ function standardized_residuals(obj::LinearModel) return r ./(std(r) .* sqrt.(1 .- h)) end +function lmplot end + function cooksleverageplot end function cooksleverageplot! end function scalelocationplot end From 6e0c603fdabd8b3d5c8d38f282e4ef6afecd5e76 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 31 Dec 2024 18:34:08 +0100 Subject: [PATCH 05/13] Minor adjustments --- ext/MakieExt.jl | 2 +- ext/StatsPlotsExt.jl | 151 ++++++++++++++++++++++++++++--------------- 2 files changed, 100 insertions(+), 53 deletions(-) diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl index 0b2cc129..bebab880 100644 --- a/ext/MakieExt.jl +++ b/ext/MakieExt.jl @@ -169,7 +169,7 @@ module MakieExt function Makie.plot!(qqp::QQPlot{<:Tuple{LinearModel}}) r = standardized_residuals(qqp[1][]) qqplot!(qqp, Normal, r, - qqline = :identity, + qqline = :fitrobust, linestyle = :dash, linewidth = 1 ) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 791f946c..0bde2a89 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -10,15 +10,19 @@ module StatsPlotsExt import GLM: scalelocationplot, scalelocationplot! import GLM: residualplot, residualplot! import GLM: residualsleverageplot, residualsleverageplot! - import GLM: quantilequantileplot, quantilequantileplot! - import StatsPlots: QQPlot, QQNorm import StatsPlots: qqplot, qqplot!, qqnorm, qqnorm! - - - - - @recipe function f(l::LinearModel) - + import GLM: lmplot + + + function lmplot(obj::LinearModel; kw...) + return plot( + residualplot(obj), + qqplot(obj), + scalelocationplot(obj), + residualsleverageplot(obj); + layout = (2,2), + kw... + ) end @userplot struct ResidualPlot{T<:Tuple{LinearModel}} @@ -70,60 +74,89 @@ module StatsPlotsExt end - @userplot struct QuantileQuantilePlot{T<:Tuple{LinearModel}} - args::T - end - #= - @recipe function f(qqp::QuantileQuantilePlot) - xlabel --> "Theoretical Quantiles" - ylabel --> "Standardized residuals" - title --> "Q-Q Residuals" - label --> "" - r = residuals(qqp.args[1]) - @series begin - seriestype := :qqnorm - linestyle := :dash - linecolor := :black - linewidth := 0.5 - r - end + function qqplot(l::LinearModel, D=Normal; + xlabel = "Theoretical Quantiles", + ylabel = "Standardized Residuals", + title = "Q-Q Residuals", + label = "", + linestyle = :dash, + linecolor = :black, + linewidth = 0.5, + kw... + ) + qqplot(D, standardized_residuals(l); + linestyle=linestyle, + linecolor=linecolor, + linewidth=linewidth, + xlabel=xlabel, + ylabel=ylabel, + title=title, + label=label, + kw...) end - =# - - @recipe function f(qqp::QuantileQuantilePlot) - xlabel --> "Theoretical Quantiles" - ylabel --> "Standardized Residuals" - title --> "Q-Q Residuals" - label --> "" - linestyle --> :dash - linecolor --> :gray - linewidth --> 0.5 - - QQPlot((Normal, standardized_residuals(qqp.args[1]))) + function qqplot!(p, l::LinearModel, D=Normal; + xlabel = "Theoretical Quantiles", + ylabel = "Standardized Residuals", + title = "Q-Q Residuals", + label = "", + linestyle = :dash, + linecolor = :black, + linewidth = 0.5, + kw... + ) + qqplot!(p, D, standardized_residuals(l); + linestyle=linestyle, + linecolor=linecolor, + linewidth=linewidth, + xlabel=xlabel, + ylabel=ylabel, + title=title, + label=label, + kw...) end - StatsPlots.recipetype(::Val{:qqplot}, obj::LinearModel, args...) = QuantileQuantilePlot((obj, args...)) - - #= - # This feels hacky but it works - qqplot(l::LinearModel, D=Normal; + function qqnorm(l::LinearModel; xlabel = "Theoretical Quantiles", ylabel = "Standardized Residuals", title = "Q-Q Residuals", label = "", linestyle = :dash, - linecolor = :gray, + linecolor = :black, linewidth = 0.5, kw... - ) = qqplot(D, standardized_residuals(l); lsty=linestyle, lcol=linecolor, lw=linewidth, xlabel=xlabel,ylabel=ylabel,title=title,label=label, kw...) - - qqplot!(p, l::LinearModel, D=Normal; kw...) = qqplot!(p, D, standardized_residuals(l); kw...) + ) + qqnorm(standardized_residuals(l); + linestyle=linestyle, + linecolor=linecolor, + linewidth=linewidth, + xlabel=xlabel, + ylabel=ylabel, + title=title, + label=label, + kw...) + end - qqnorm(l::LinearModel; kw...) = qqnorm(standardized_residuals(l); kw...) - qqnorm!(p, l::LinearModel; kw...) = qqnorm!(p, standardized_residuals(l), kw...) - =# + function qqnorm!(p, l::LinearModel; + xlabel = "Theoretical Quantiles", + ylabel = "Standardized Residuals", + title = "Q-Q Residuals", + label = "", + linestyle = :dash, + linecolor = :black, + linewidth = 0.5, + kw...) + qqnorm!(p, standardized_residuals(l), + linestyle=linestyle, + linecolor=linecolor, + linewidth=linewidth, + xlabel=xlabel, + ylabel=ylabel, + title=title, + label=label, + kw...) + end @userplot struct ResidualsLeveragePlot{T<:Tuple{LinearModel}} args::T @@ -141,6 +174,7 @@ module StatsPlotsExt ymax = maximum(abs.(r)) ylims --> (-1.1*ymax, 1.1*ymax) + xlims --> (0.0, :auto) @series begin seriestype := :scatter @@ -165,7 +199,7 @@ module StatsPlotsExt if !isempty(cook_levels) cookfun = (h,D,k) -> sqrt(D * k * (1-h) / h) xmin, xmax = extrema(h) - xs = LinRange(xmin, 1.1*xmax, 50) + xs = LinRange(xmin, 1.2*xmax, 50) for D in cook_levels @series begin seriestype := :path @@ -178,17 +212,30 @@ module StatsPlotsExt @series begin seriestype := :scatter markeralpha := 0.0 - series_annotation := [("$D", 9, :gray)] + series_annotation := [("$D", 8, :gray)] [1.02 * xs[end]], [cookfun(xs[end],D,k)] end + @series begin + seriestype := :scatter + markeralpha := 0.0 + series_annotation := [("$D", 8, :gray)] + [1.02 * xs[end]], [-cookfun(xs[end],D,k)] + end @series begin seriestype := :path linecolor := :gray linestyle := :dash - xs, -cookfun.(xs,D,k) + xs, -cookfun.(xs,D,k) end end + begin + seriestype := :path + linecolor := :gray + linestyle := :dash + label := "Cook's Distance" + [-1,-1],[-1,-1] + end end From df5a618c1c08fa1ab3c5d6939863a4ae6dc05863 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Wed, 1 Jan 2025 22:15:07 +0100 Subject: [PATCH 06/13] Add docstrings --- ext/StatsPlotsExt.jl | 83 +++++++++++++++++++++++++++----------------- src/plots.jl | 70 +++++++++++++++++++++++++++++++++++-- 2 files changed, 119 insertions(+), 34 deletions(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 0bde2a89..944bc2b2 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -29,7 +29,12 @@ module StatsPlotsExt args::T end - @recipe function f(rp::ResidualPlot; ) + @recipe function f(rp::ResidualPlot; + axislines = true, + axislinestyle = :dot, + axislinecolor = :black, + axislinewidth = 0.5 + ) xlabel --> "Fitted values" ylabel --> "Residuals" title --> "Residuals vs Fitted" @@ -42,14 +47,15 @@ module StatsPlotsExt seriestype := :scatter y, r end - - @series begin - seriestype := :hline - linecolor := :black - linestyle := :dash - linewidth := 0.5 - label := "" - [0.0] + if axislines + @series begin + seriestype := :hline + linecolor := axislinecolor + linestyle := axislinestyle + linewidth := axislinewidth + label := "" + [0.0] + end end nothing end @@ -162,7 +168,16 @@ module StatsPlotsExt args::T end - @recipe function f(rlp::ResidualsLeveragePlot; cook_levels = [0.5, 1.0]) + @recipe function f(rlp::ResidualsLeveragePlot; + axislines = true, + axislinestyle = :dot, + axislinecolor = :gray, + axislinewidth = 1, + cookslevels = [0.5, 1.0], + cookslinecolor = :gray, + cookslinestyle = :dash, + cookslinewidth = 1 + ) xlabel --> "Leverage" ylabel --> "Standardized residuals" title --> "Residuals vs Leverage" @@ -180,31 +195,33 @@ module StatsPlotsExt seriestype := :scatter h, r end - - @series begin - seriestype := :hline - linecolor := :black - linestyle := :dash - linewidth := 0.5 - label := "" - [0.0] - end - @series begin - seriestype := :vline - linecolor := :black - linestyle := :dash - linewidth := 0.5 - [0.0] + if axislines + @series begin + seriestype := :hline + linecolor := axislinecolor + linestyle := axislinestyle + linewidth := axislinewidth + label := "" + [0.0] + end + @series begin + seriestype := :vline + linecolor := axislinecolor + linestyle := axislinestyle + linewidth := axislinewidth + [0.0] + end end - if !isempty(cook_levels) + if !isempty(cookslevels) cookfun = (h,D,k) -> sqrt(D * k * (1-h) / h) xmin, xmax = extrema(h) xs = LinRange(xmin, 1.2*xmax, 50) for D in cook_levels @series begin seriestype := :path - linecolor := :gray - linestyle := :dash + linecolor := cookslinecolor + linestyle := cookslinestyle + linewidth := cookslinewidth annotations := [(xs[end], cookfun(xs[end],D,k), ("$D", 8, :gray))] xs, cookfun.(xs,D,k) end @@ -224,15 +241,17 @@ module StatsPlotsExt @series begin seriestype := :path - linecolor := :gray - linestyle := :dash + linecolor := cookslinecolor + linestyle := cookslinestyle + linewidth := cookslinewidth xs, -cookfun.(xs,D,k) end end begin seriestype := :path - linecolor := :gray - linestyle := :dash + linecolor := cookslinecolor + linestyle := cookslinestyle + linewidth := cookslinewidth label := "Cook's Distance" [-1,-1],[-1,-1] end diff --git a/src/plots.jl b/src/plots.jl index 29eece21..0f176e8e 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -14,15 +14,81 @@ function standardized_residuals(obj::LinearModel) return r ./(std(r) .* sqrt.(1 .- h)) end + +""" + lmplot(obj::LinearModel; kw...) + +Display several summary plots of a linear model. + +Keyword arguments for the plotting backend such as `size` are supported. If using Makie, only keyword arguments to `Figure` are supported. + +## Examples +```julia-repl +julia> using GLM, StatsPlots + +julia> X = randn(30, 5); y = X * randn(5) + 0.3*randn(30) + +julia> l = lm(X,y) + +julia> lmplot(l) +``` +""" function lmplot end +""" + cooksleverageplot(obj::LinearModel; kw...) + +Plot the Cook's distances of a linear model against its leverages. + +Keyword arguments are passed to the plotting backend. +""" function cooksleverageplot end function cooksleverageplot! end + +""" + scalelocationplot(obj::LinearModel, kw...) + +Plot the root standardized residuals of a linear model against its fitted values. + +Keyword arguments are passed to the plotting backend. +""" function scalelocationplot end function scalelocationplot! end + +""" + residualplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its fitted values. + +## keyword arguments + +* `axislines = true` whether to display a line on the x axis. +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` + +Other keyword arguments are passed to the plotting backend. +""" function residualplot end function residualplot! end + +""" + residualsleverageplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its leverages. + +## keyword arguments + +* `axislines = true` +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` +* `cookslevels = [0.5,2.0]` Levels curves of Cook's distance to display. +* `cookslinecolor` +* `cookslinestyle` +* `cookslinewidth` + +Other keyword arguments are passed to the plotting backend. +""" function residualsleverageplot end function residualsleverageplot! end -function quantilequantileplot end -function quantilequantileplot! end From 45c3ff213b26151485db72f45f295c5ee949d377 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Wed, 1 Jan 2025 23:59:46 +0100 Subject: [PATCH 07/13] Fix typo --- ext/StatsPlotsExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 944bc2b2..5c1794d6 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -216,7 +216,7 @@ module StatsPlotsExt cookfun = (h,D,k) -> sqrt(D * k * (1-h) / h) xmin, xmax = extrema(h) xs = LinRange(xmin, 1.2*xmax, 50) - for D in cook_levels + for D in cookslevels @series begin seriestype := :path linecolor := cookslinecolor From 1846e2352f8ab8ad0db14581f1b3e3182bee4dc5 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 7 Jan 2025 09:57:21 +0100 Subject: [PATCH 08/13] Add backwards compatibility --- Project.toml | 6 ++++-- src/GLM.jl | 11 +++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index d3d93c31..a8f843cb 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -18,12 +19,12 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [weakdeps] -StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" [extensions] -StatsPlotsExt = "StatsPlots" MakieExt = "Makie" +StatsPlotsExt = "StatsPlots" [compat] CSV = "0.7, 0.8, 0.9, 0.10" @@ -32,6 +33,7 @@ DataFrames = "0.22, 1" Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" RDatasets = "0.5, 0.6, 0.7" Reexport = "0.1, 0.2, 1.0" +Requires = "1.3.0" SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1, 2.0" Statistics = "1" StatsAPI = "1.4" diff --git a/src/GLM.jl b/src/GLM.jl index 72502dae..cdc02f4f 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -106,6 +106,17 @@ module GLM pivoted_qr!(A; kwargs...) = qr!(A, ColumnNorm(); kwargs...) end + if !isdefined(Base, :get_extension) + using Requires + end + + function __init__() + @static if !isdefined(Base, :get_extension) + @require StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" include("../ext/StatsPlotsExt.jl") + @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("../ext/MakieExt.jl") + end + end + const COMMON_FIT_KWARGS_DOCS = """ - `dropcollinear::Bool=true`: Controls whether or not a model matrix less-than-full rank is accepted. From b70a682b8c5bd5e185185e04751b6630984b5eae Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 7 Jan 2025 12:04:08 +0100 Subject: [PATCH 09/13] Add compat --- Project.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index a8f843cb..3242abf7 100644 --- a/Project.toml +++ b/Project.toml @@ -31,6 +31,7 @@ CSV = "0.7, 0.8, 0.9, 0.10" CategoricalArrays = "0.8, 0.9, 0.10" DataFrames = "0.22, 1" Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25" +Makie = "0.21" RDatasets = "0.5, 0.6, 0.7" Reexport = "0.1, 0.2, 1.0" Requires = "1.3.0" @@ -49,8 +50,10 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] From fd080fe99a5e21352d498ee4f5e7d17b45409ab8 Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 25 Feb 2025 22:22:50 +0100 Subject: [PATCH 10/13] Add simple tests --- test/plots.jl | 52 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 4 ++++ 2 files changed, 56 insertions(+) create mode 100644 test/plots.jl diff --git a/test/plots.jl b/test/plots.jl new file mode 100644 index 00000000..5736e9a1 --- /dev/null +++ b/test/plots.jl @@ -0,0 +1,52 @@ +using Test +using GLM +using StatsPlots +using GLM: standardized_residuals, leverage + + +@testset "Utility functions" begin + rng = StableRNG(2025) + X = randn(rng,10,3) + y = randn(rng,10) + l = lm(X,y) + h = leverage(l) + r = standardized_residuals(l) + @test all(!isnan, r) + @test all(>=(0.0), h) +end + +@testset "StatsPlots Recipes" begin + # NB. These tests follow the tests of StatsPlots. They mostly check that the functions don't crash + rng = StableRNG(2025) + X = randn(rng,10,3) + y = randn(rng,10) + l = lm(X,y) + @testset "residualplot" begin + pl = residualplot(l) + @test show(devnull, pl) isa Nothing + end + @testset "residualsleverageplot" begin + pl = residualsleverageplot(l) + @test show(devnull, pl) isa Nothing + end + @testset "scalelocationplot" begin + pl = scalelocationplot(l) + @test show(devnull, pl) isa Nothing + end + @testset "qqplot" begin + pl = qqplot(l) + @test show(devnull, pl) isa Nothing + end + @testset "cooksleverageplot" begin + pl = cooksleverageplot(l) + @test show(devnull, pl) isa Nothing + end + @testset "lmplot" begin + pl = lmplot(l) + @test show(devnull, pl) isa Nothing + end +end + +@testset "Makie Recipes" begin + +end diff --git a/test/runtests.jl b/test/runtests.jl index fb0115f5..32bae57f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2044,3 +2044,7 @@ end @test coef(ft) ≈ [9.648767705301294, -0.11274823562143056, 0.1907889126252095, -0.8123086879222496] @test_throws DomainError glm(@formula(Column1 ~ Column2 + Column3 + Column4), df, Gamma(), LogLink(), start = fill(NaN, 4)) end + +@testset "Plots" begin + include("plots.jl") +end From 3e3e294470cf6b531725ac56d2047c4de0fc2f0d Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Tue, 25 Feb 2025 22:32:44 +0100 Subject: [PATCH 11/13] Add StatsPlots as a test dependency --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3242abf7..541c41ac 100644 --- a/Project.toml +++ b/Project.toml @@ -57,4 +57,4 @@ StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["CategoricalArrays", "CSV", "DataFrames", "Downloads", "RDatasets", "StableRNGs", "Test"] +test = ["CategoricalArrays", "CSV", "DataFrames", "Downloads", "RDatasets", "StableRNGs", "StatsPlots", "Test"] From 4d060dc3869a8a83fc9338af091741804efb43da Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Sat, 1 Mar 2025 22:15:11 +0100 Subject: [PATCH 12/13] Add submodules for plot recipes --- ext/MakieExt.jl | 10 ++--- ext/StatsPlotsExt.jl | 14 +++---- src/GLM.jl | 8 +--- src/makie.jl | 89 ++++++++++++++++++++++++++++++++++++++++++++ src/plots.jl | 79 --------------------------------------- src/statsplots.jl | 89 ++++++++++++++++++++++++++++++++++++++++++++ test/plots.jl | 63 ++++++++++++++++++++++++++++--- 7 files changed, 249 insertions(+), 103 deletions(-) create mode 100644 src/makie.jl create mode 100644 src/statsplots.jl diff --git a/ext/MakieExt.jl b/ext/MakieExt.jl index bebab880..8266d3ed 100644 --- a/ext/MakieExt.jl +++ b/ext/MakieExt.jl @@ -3,11 +3,11 @@ module MakieExt using GLM using Distributions using GLM: leverage, standardized_residuals - import GLM: cooksleverageplot, cooksleverageplot! - import GLM: scalelocationplot, scalelocationplot! - import GLM: residualplot, residualplot! - import GLM: residualsleverageplot, residualsleverageplot! - import GLM: lmplot + import GLM.MakieRecipes: cooksleverageplot, cooksleverageplot! + import GLM.MakieRecipes: scalelocationplot, scalelocationplot! + import GLM.MakieRecipes: residualplot, residualplot! + import GLM.MakieRecipes: residualsleverageplot, residualsleverageplot! + import GLM.MakieRecipes: lmplot import Makie: qqplot, qqplot! diff --git a/ext/StatsPlotsExt.jl b/ext/StatsPlotsExt.jl index 5c1794d6..31df2106 100644 --- a/ext/StatsPlotsExt.jl +++ b/ext/StatsPlotsExt.jl @@ -6,12 +6,12 @@ module StatsPlotsExt using Distributions using GLM: leverage, standardized_residuals using RecipesBase: recipetype - import GLM: cooksleverageplot, cooksleverageplot! - import GLM: scalelocationplot, scalelocationplot! - import GLM: residualplot, residualplot! - import GLM: residualsleverageplot, residualsleverageplot! + import GLM.PlotsRecipes: cooksleverageplot, cooksleverageplot! + import GLM.PlotsRecipes: scalelocationplot, scalelocationplot! + import GLM.PlotsRecipes: residualplot, residualplot! + import GLM.PlotsRecipes: residualsleverageplot, residualsleverageplot! + import GLM.PlotsRecipes: lmplot import StatsPlots: qqplot, qqplot!, qqnorm, qqnorm! - import GLM: lmplot function lmplot(obj::LinearModel; kw...) @@ -102,7 +102,7 @@ module StatsPlotsExt kw...) end - function qqplot!(p, l::LinearModel, D=Normal; + function qqplot!(p::RecipesBase.AbstractPlot, l::LinearModel, D=Normal; xlabel = "Theoretical Quantiles", ylabel = "Standardized Residuals", title = "Q-Q Residuals", @@ -144,7 +144,7 @@ module StatsPlotsExt kw...) end - function qqnorm!(p, l::LinearModel; + function qqnorm!(p::RecipesBase.AbstractPlot, l::LinearModel; xlabel = "Theoretical Quantiles", ylabel = "Standardized Residuals", title = "Q-Q Residuals", diff --git a/src/GLM.jl b/src/GLM.jl index cdc02f4f..6b338a19 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -67,12 +67,6 @@ module GLM predict, # make predictions ftest # compare models with an F test - # Plot functions - export cooksleverageplot, cooksleverageplot! - export scalelocationplot, scalelocationplot! - export residualplot, residualplot! - export residualsleverageplot, residualsleverageplot! - export lmplot const FP = AbstractFloat @@ -152,5 +146,7 @@ module GLM include("negbinfit.jl") include("deprecated.jl") include("plots.jl") + include("statsplots.jl") + include("makie.jl") end # module diff --git a/src/makie.jl b/src/makie.jl new file mode 100644 index 00000000..5a89837c --- /dev/null +++ b/src/makie.jl @@ -0,0 +1,89 @@ +module MakieRecipes + + # Plot functions + export cooksleverageplot, cooksleverageplot! + export scalelocationplot, scalelocationplot! + export residualplot, residualplot! + export residualsleverageplot, residualsleverageplot! + export lmplot + + +""" + lmplot(obj::LinearModel; kw...) + +Display several summary plots of a linear model. + +Keyword arguments for the plotting backend such as `size` are supported. If using Makie, only keyword arguments to `Figure` are supported. + +## Examples +```julia-repl +julia> using GLM, Makie, GLM.MakieRecipes + +julia> X = randn(30, 5); y = X * randn(5) + 0.3*randn(30) + +julia> l = lm(X,y) + +julia> lmplot(l) +``` +""" +function lmplot end + +""" + cooksleverageplot(obj::LinearModel; kw...) + +Plot the Cook's distances of a linear model against its leverages. + +Keyword arguments are passed to the plotting backend. +""" +function cooksleverageplot end +function cooksleverageplot! end + +""" + scalelocationplot(obj::LinearModel, kw...) + +Plot the root standardized residuals of a linear model against its fitted values. + +Keyword arguments are passed to the plotting backend. +""" +function scalelocationplot end +function scalelocationplot! end + +""" + residualplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its fitted values. + +## keyword arguments + +* `axislines = true` whether to display a line on the x axis. +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` + +Other keyword arguments are passed to the plotting backend. +""" +function residualplot end +function residualplot! end + +""" + residualsleverageplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its leverages. + +## keyword arguments + +* `axislines = true` +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` +* `cookslevels = [0.5,2.0]` Levels curves of Cook's distance to display. +* `cookslinecolor` +* `cookslinestyle` +* `cookslinewidth` + +Other keyword arguments are passed to the plotting backend. +""" +function residualsleverageplot end +function residualsleverageplot! end + +end diff --git a/src/plots.jl b/src/plots.jl index 0f176e8e..469ce548 100644 --- a/src/plots.jl +++ b/src/plots.jl @@ -13,82 +13,3 @@ function standardized_residuals(obj::LinearModel) h = leverage(obj) return r ./(std(r) .* sqrt.(1 .- h)) end - - -""" - lmplot(obj::LinearModel; kw...) - -Display several summary plots of a linear model. - -Keyword arguments for the plotting backend such as `size` are supported. If using Makie, only keyword arguments to `Figure` are supported. - -## Examples -```julia-repl -julia> using GLM, StatsPlots - -julia> X = randn(30, 5); y = X * randn(5) + 0.3*randn(30) - -julia> l = lm(X,y) - -julia> lmplot(l) -``` -""" -function lmplot end - -""" - cooksleverageplot(obj::LinearModel; kw...) - -Plot the Cook's distances of a linear model against its leverages. - -Keyword arguments are passed to the plotting backend. -""" -function cooksleverageplot end -function cooksleverageplot! end - -""" - scalelocationplot(obj::LinearModel, kw...) - -Plot the root standardized residuals of a linear model against its fitted values. - -Keyword arguments are passed to the plotting backend. -""" -function scalelocationplot end -function scalelocationplot! end - -""" - residualplot(obj::LinearModel, kw...) - -Plot the residuals of a linear model against its fitted values. - -## keyword arguments - -* `axislines = true` whether to display a line on the x axis. -* `axislinecolor` -* `axislinestyle` -* `axislinewidth` - -Other keyword arguments are passed to the plotting backend. -""" -function residualplot end -function residualplot! end - -""" - residualsleverageplot(obj::LinearModel, kw...) - -Plot the residuals of a linear model against its leverages. - -## keyword arguments - -* `axislines = true` -* `axislinecolor` -* `axislinestyle` -* `axislinewidth` -* `cookslevels = [0.5,2.0]` Levels curves of Cook's distance to display. -* `cookslinecolor` -* `cookslinestyle` -* `cookslinewidth` - -Other keyword arguments are passed to the plotting backend. -""" -function residualsleverageplot end -function residualsleverageplot! end diff --git a/src/statsplots.jl b/src/statsplots.jl new file mode 100644 index 00000000..b9147966 --- /dev/null +++ b/src/statsplots.jl @@ -0,0 +1,89 @@ +module PlotsRecipes + + # Plot functions + export cooksleverageplot, cooksleverageplot! + export scalelocationplot, scalelocationplot! + export residualplot, residualplot! + export residualsleverageplot, residualsleverageplot! + export lmplot + + +""" + lmplot(obj::LinearModel; kw...) + +Display several summary plots of a linear model. + +Keyword arguments for the plotting backend such as `size` are supported. If using Makie, only keyword arguments to `Figure` are supported. + +## Examples +```julia-repl +julia> using GLM, StatsPlots, GLM.PlotsRecipes + +julia> X = randn(30, 5); y = X * randn(5) + 0.3*randn(30) + +julia> l = lm(X,y) + +julia> lmplot(l) +``` +""" +function lmplot end + +""" + cooksleverageplot(obj::LinearModel; kw...) + +Plot the Cook's distances of a linear model against its leverages. + +Keyword arguments are passed to the plotting backend. +""" +function cooksleverageplot end +function cooksleverageplot! end + +""" + scalelocationplot(obj::LinearModel, kw...) + +Plot the root standardized residuals of a linear model against its fitted values. + +Keyword arguments are passed to the plotting backend. +""" +function scalelocationplot end +function scalelocationplot! end + +""" + residualplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its fitted values. + +## keyword arguments + +* `axislines = true` whether to display a line on the x axis. +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` + +Other keyword arguments are passed to the plotting backend. +""" +function residualplot end +function residualplot! end + +""" + residualsleverageplot(obj::LinearModel, kw...) + +Plot the residuals of a linear model against its leverages. + +## keyword arguments + +* `axislines = true` +* `axislinecolor` +* `axislinestyle` +* `axislinewidth` +* `cookslevels = [0.5,2.0]` Levels curves of Cook's distance to display. +* `cookslinecolor` +* `cookslinestyle` +* `cookslinewidth` + +Other keyword arguments are passed to the plotting backend. +""" +function residualsleverageplot end +function residualsleverageplot! end + +end diff --git a/test/plots.jl b/test/plots.jl index 5736e9a1..b14a63d4 100644 --- a/test/plots.jl +++ b/test/plots.jl @@ -1,7 +1,11 @@ using Test using GLM using StatsPlots +using CairoMakie +#using Makie using GLM: standardized_residuals, leverage +import GLM.PlotsRecipes +import GLM.MakieRecipes @testset "Utility functions" begin @@ -22,31 +26,78 @@ end y = randn(rng,10) l = lm(X,y) @testset "residualplot" begin - pl = residualplot(l) + pl = PlotsRecipes.residualplot(l) + @test show(devnull, pl) isa Nothing + PlotsRecipes.residualplot!(pl, l) @test show(devnull, pl) isa Nothing end @testset "residualsleverageplot" begin - pl = residualsleverageplot(l) + pl = PlotsRecipes.residualsleverageplot(l) + @test show(devnull, pl) isa Nothing + PlotsRecipes.residualsleverageplot!(pl, l) @test show(devnull, pl) isa Nothing end @testset "scalelocationplot" begin - pl = scalelocationplot(l) + pl = PlotsRecipes.scalelocationplot(l) + @test show(devnull, pl) isa Nothing + PlotsRecipes.scalelocationplot!(pl, l) @test show(devnull, pl) isa Nothing end @testset "qqplot" begin - pl = qqplot(l) + pl = StatsPlots.qqplot(l) + @test show(devnull, pl) isa Nothing + StatsPlots.qqplot!(pl, l) + @test show(devnull, pl) isa Nothing + end + @testset "qqnorm" begin + pl = StatsPlots.qqnorm(l) + @test show(devnull, pl) isa Nothing + StatsPlots.qqnorm!(pl, l) @test show(devnull, pl) isa Nothing end + @testset "cooksleverageplot" begin - pl = cooksleverageplot(l) + pl = PlotsRecipes.cooksleverageplot(l) + @test show(devnull, pl) isa Nothing + PlotsRecipes.cooksleverageplot!(pl, l) @test show(devnull, pl) isa Nothing end @testset "lmplot" begin - pl = lmplot(l) + pl = PlotsRecipes.lmplot(l) @test show(devnull, pl) isa Nothing end end @testset "Makie Recipes" begin + rng = StableRNG(2025) + X = randn(rng, 10, 3) + y = randn(rng, 10) + l = lm(X,y) + @testset "residualplot" begin + fig, ax, plt = MakieRecipes.residualplot(l) + @test plt isa Makie.Plot + end + @testset "residualsleverageplot" begin + fig, ax, plt = MakieRecipes.residualsleverageplot(l) + @test plt isa Makie.Plot + end + @testset "scalelocationplot" begin + fig, ax, plt = MakieRecipes.scalelocationplot(l) + @test plt isa Makie.Plot + end + @testset "qqplot" begin + fig, ax, plt = CairoMakie.qqplot(l) + @test plt isa Makie.Plot + end + @testset "cooksleverageplot" begin + fig, ax, plt = MakieRecipes.cooksleverageplot(l) + @test plt isa Makie.Plot + end + @testset "lmplot" begin + fig = MakieRecipes.lmplot(l) + @test fig isa Makie.Figure + + end + end From ae37c1cbbc03087b18e378c845aa57785aaedadb Mon Sep 17 00:00:00 2001 From: irregular-rhomboid Date: Sat, 1 Mar 2025 22:25:21 +0100 Subject: [PATCH 13/13] Update Project.toml --- Project.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 541c41ac..98a30861 100644 --- a/Project.toml +++ b/Project.toml @@ -48,6 +48,7 @@ julia = "1.6" [extras] CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @@ -57,4 +58,4 @@ StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["CategoricalArrays", "CSV", "DataFrames", "Downloads", "RDatasets", "StableRNGs", "StatsPlots", "Test"] +test = ["CategoricalArrays", "CSV", "CairoMakie", "DataFrames", "Downloads", "RDatasets", "StableRNGs", "StatsPlots", "Test"]