VIF

Project.toml

@@ -1,10 +1,11 @@
 name = "StatsAPI"
 uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 authors = ["Milan Bouchet-Valat <[email protected]"]
-version = "1.6.0"
+version = "1.7.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 julia = "1"

src/StatsAPI.jl

@@ -1,6 +1,7 @@
 module StatsAPI
 
 using LinearAlgebra
+using Statistics
 
 include("statisticalmodel.jl")
 include("regressionmodel.jl")

src/regressionmodel.jl

@@ -132,3 +132,43 @@ function linearpredictor end
 In-place version of [`linearpredictor`](@ref), storing the result in `storage`.
 """
 function linearpredictor! end
+
+"""
+    vif(m::RegressionModel)
+
+Compute the variance inflation factor (VIF).
+
+Returns a vector of inflation factors computed for each coefficient except
+for the intercept.
+In other words, the corresponding coefficients are `coefnames(m)[2:end]`.
+
+The [variance inflation factor (VIF)](https://en.wikipedia.org/wiki/Variance_inflation_factor) measures
+the increase in the variance of a parameter's estimate in a model with multiple parameters relative to
+the variance of a parameter's estimate in a model containing only that parameter.
+
+See also [`coefnames`](@ref), [`gvif`](@ref).
+
+!!! warning
+    This method will fail if there is (numerically) perfect multicollinearity,
+    i.e. rank deficiency. In that case though, the VIF
+    is not particularly informative anyway.
+"""
+function vif(model::RegressionModel)
+    mm = Statistics.cov2cor!(vcov(model), stderror(model))
+    # TODO: replace with hasintercept() when implemented (xref #17)
+    all(==(1), view(modelmatrix(model), :, 1)) ||
+        throw(ArgumentError("VIF only defined for models with an intercept term"))
+    mm = @view mm[2:end, 2:end]
+    size(mm, 2) > 1 ||
+        throw(ArgumentError("VIF not meaningful for models with only one non-intercept term"))
+    # NB: The correlation matrix is positive definite and hence invertible
+    #     unless there is perfect rank deficiency, hence the warning.
+    # so we want diag(inv(mm)) but directly computing inverses is bad.
+    # that said, these are typically small-ish matrices and this is Simple.
+    return diag(inv(mm))
+end
+
+# This generic function is owned by StatsModels.jl, which is the sole provider
+# of the default definition. StatsModels needs to own the default definition
+# because it depends on the @formula interface.
+function gvif end

test/regressionmodel.jl

@@ -1,18 +1,43 @@
 module TestRegressionModel
 
 using Test, LinearAlgebra, StatsAPI
-using StatsAPI: RegressionModel, crossmodelmatrix
+using StatsAPI: RegressionModel, crossmodelmatrix, vif
 
 struct MyRegressionModel <: RegressionModel
 end
 
 StatsAPI.modelmatrix(::MyRegressionModel) = [1 2; 3 4]
+StatsAPI.vcov(::MyRegressionModel) = [1 0; 0 1]
+
+struct MyRegressionModel2 <: RegressionModel
+end
+
+StatsAPI.modelmatrix(::MyRegressionModel2) = [1 2; 1 2]
+StatsAPI.vcov(::MyRegressionModel2) = [1 0; 0 1]
+
+struct MyRegressionModel3 <: RegressionModel
+end
+
+StatsAPI.modelmatrix(::MyRegressionModel3) = [1 2 3; 1 2 3]
+StatsAPI.vcov(::MyRegressionModel3) = [1 0 0; 0 1 0; 0 0 1]
+
 
 @testset "TestRegressionModel" begin
     m = MyRegressionModel()
 
     @test crossmodelmatrix(m) == [10 14; 14 20]
     @test crossmodelmatrix(m) isa Symmetric
+
+    # no intercept term
+    @test_throws ArgumentError vif(m)
+
+    m2 = MyRegressionModel2()
+    # only one non intercept term
+    @test_throws ArgumentError vif(m2)
+
+    m3 = MyRegressionModel3()
+    # vcov is identity, so the VIF is just 1
+    @test vif(m3) ≈ [1, 1]
 end
 
-end # module TestRegressionModel
+end # module TestRegressionModel

test/statisticalmodel.jl

@@ -36,4 +36,4 @@ StatsAPI.nobs(::MyStatisticalModel) = 100
     @test adjr2 === adjr²
 end
 
-end # module TestStatisticalModel
+end # module TestStatisticalModel