Try w/ “real” data:
using DataFrames, GLM, Statistics, LinearAlgebra, RDatasets
airquality = rename(dataset("datasets", "airquality"), "Solar.R" => "Solar_R")
#
y = airquality.Wind
X1 = airquality.Temp
X2 = airquality.Solar_R
X3 = airquality.Ozone
d = DataFrame(X1=X1, X2=X2, X3=X3, y=y)
d = d[completecases(d), :]
X = [one.(d.X1) d.X1 d.X2 d.X3] .|> Float64
y = d.y .|> Float64
#
β_LS = X \ y
vifm = diag(inv(cor(X[:,2:end])))
#
m_y = lm(@formula(y ~ 1 + X1 + X2 + X3), d)
m_X1 = lm(@formula(X1 ~ 1 + X2 + X3), d)
m_X2 = lm(@formula(X2 ~ 1 + X1 + X3), d)
m_X3 = lm(@formula(X3 ~ 1 + X1 + X2), d)
vif_X1 = 1.0/(1.0-r2(m_X1))
vif_X2 = 1.0/(1.0-r2(m_X2))
vif_X3 = 1.0/(1.0-r2(m_X3))
#
vifm ≈ [vif_X1; vif_X2; vif_X3;]
#
vif_GLM(mod) = diag(inv(cor(mod.model.pp.X[:,2:end])))
vif_GLM(m_y) ≈ vifm ≈ [vif_X1; vif_X2; vif_X3;]