Cooks Distance tests

Project.toml

@@ -23,9 +23,9 @@ Distributions = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24"
 RDatasets = "0.5, 0.6, 0.7"
 Reexport = "0.1, 0.2, 1.0"
 SpecialFunctions = "0.6, 0.7, 0.8, 0.9, 0.10, 1"
-StatsBase = "0.30, 0.31, 0.32, 0.33"
+StatsBase = "0.33.5"
 StatsFuns = "0.6, 0.7, 0.8, 0.9"
-StatsModels = "0.6"
+StatsModels = "0.6.22"
 julia = "1"
 
 [extras]

src/GLM.jl

@@ -17,7 +17,8 @@ module GLM
     import SpecialFunctions: erfc, erfcinv, digamma, trigamma
     export coef, coeftable, confint, deviance, nulldeviance, dof, dof_residual,
            loglikelihood, nullloglikelihood, nobs, stderror, vcov, residuals, predict,
-           fitted, fit, fit!, model_response, response, modelmatrix, r2, r², adjr2, adjr²
+           fitted, fit, fit!, model_response, response, modelmatrix, r2, r², adjr2, adjr², 
+           cooksdistance
 
     export
         # types

src/lm.jl

@@ -243,3 +243,29 @@ function confint(obj::LinearModel; level::Real=0.95)
     hcat(coef(obj),coef(obj)) + stderror(obj) *
     quantile(TDist(dof_residual(obj)), (1. - level)/2.) * [1. -1.]
 end
+
+"""
+    cooksdistance(obj::LinearModel)
+
+Compute [Cook's distance](https://en.wikipedia.org/wiki/Cook%27s_distance)
+for each observation in linear model `obj`, giving an estimate of the influence
+of each data point.
+Currently only implemented for LinearModel models without weights.
+"""
+function StatsBase.cooksdistance(obj::LinearModel)
+    u = residuals(obj)
+    mse = dispersion(obj,true)
+    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."))
+    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
+end

test/runtests.jl

@@ -49,6 +49,38 @@ linreg(x::AbstractVecOrMat, y::AbstractVector) = qr!(simplemm(x)) \ y
     @test coef(lm3) == coef(lm4) ≈ [11, 1, 2, 3, 4]
 end
 
+@testset "Linear Model Cooks Distance - refers to PR #368 and issue #414" begin
+    st_df = DataFrame( 
+        Y=[6.4, 7.4, 10.4, 15.1, 12.3 , 11.4],
+        XA=[1.5, 6.5, 11.5, 19.9, 17.0, 15.5],
+        XB=[1.8, 7.8, 11.8, 20.5, 17.3, 15.8], 
+        XC=[3., 13., 23., 39.8, 34., 31.],
+        # values from SAS proc reg
+        CooksD_base=[1.4068501943, 0.176809102, 0.0026655177, 1.0704009915, 0.0875726457, 0.1331183932], 
+        CooksD_noint=[0.0076891801, 0.0302993877, 0.0410262965, 0.0294348488, 0.0691589296, 0.0273045538], 
+        CooksD_multi=[1.7122291956,	18.983407026, 0.000118078, 0.8470797843, 0.0715921999, 0.1105843157 ],
+        )
+
+    # linear regression
+    t_lm_base = lm(@formula(Y ~ XA), st_df)
+    @test isapprox(st_df.CooksD_base, cooksdistance(t_lm_base))
+
+    # linear regression, no intercept 
+    t_lm_noint = lm(@formula(Y ~ XA +0), st_df)
+    @test isapprox(st_df.CooksD_noint, cooksdistance(t_lm_noint))
+
+    # linear regression, two collinear variables (Variance inflation factor ≊ 250)
+    t_lm_multi = lm(@formula(Y ~ XA + XB), st_df)
+    @test isapprox(st_df.CooksD_multi, cooksdistance(t_lm_multi))
+
+    # linear regression, two full collinear variables (XC = 2 XA) hence should get the same results as the original
+    # after pivoting
+    t_lm_colli = lm(@formula(Y ~ XA + XC), st_df, dropcollinear=true)
+    ## Currently the test fails as the collinear variable is not droped from a ```modelmatrix(obj)``` call.
+    @test_throws ArgumentError isapprox(st_df.CooksD_base, cooksdistance(t_lm_colli))
+
+end
+
 @testset "linear model with weights" begin 
     df = dataset("quantreg", "engel")
     N = nrow(df)