Add Basic Plots Recipes

JuliaStats · irregular-rhomboid · Dec 8, 2024 · Dec 24, 2024 · Dec 26, 2024 · Dec 31, 2024
commit 45c5cd914fd0afb8716ea4e199f706738575cce6
diff --git 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
@@ -0,0 +1,3 @@
+module MakieExt
+
+end
diff --git a/ext/StatsPlotsExt.jl 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
@@ -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
@@ -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