Adding frequentist numerical tests + Integrate StatsAPI

Project.toml

@@ -12,6 +12,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"

src/CRRao.jl

@@ -3,6 +3,7 @@ module CRRao
 using DataFrames, GLM, Turing, StatsModels, StatsBase
 using StatsBase, Distributions, LinearAlgebra
 using Optim, NLSolversBase, Random
+import StatsBase: coeftable, r2, adjr2, loglikelihood, aic, bic, predict, residuals, cooksdistance, fit
 
 """
 ```julia

src/fitmodel.jl

@@ -5,7 +5,7 @@ FrequentistRegression{RegressionType}
 
 Type to represent frequentist regression models returned by `fit` functions. This type is used internally by the package to represent all frequentist regression models. `RegressionType` is a `Symbol` representing the model class.
 """
-struct FrequentistRegression{RegressionType}
+struct FrequentistRegression{RegressionType} <: RegressionModel
     model
     formula::FormulaTerm
     link
@@ -27,7 +27,7 @@ BayesianRegression{RegressionType}
 
 Type to represent bayesian regression models returned by `fit` functions. This type is used internally by the package to represent all bayesian regression models. `RegressionType` is a `Symbol` representing the model class.
 """
-struct BayesianRegression{RegressionType}
+struct BayesianRegression{RegressionType} <: RegressionModel
     chain
     formula::FormulaTerm
     link

src/frequentist/getter.jl

@@ -221,7 +221,7 @@ function predict(container::FrequentistRegression{:LogisticRegression}, newdata:
         return Probit_Link.(z)
     elseif (container.link == GLM.CauchitLink)
         return Cauchit_Link.(z)
-    elseif (container.link == GLM.Cloglog)
+    elseif (container.link == GLM.CloglogLink)
         return Cloglog_Link.(z)
     end
 end

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"

test/numerical/frequentist/LinearRegression.jl

@@ -0,0 +1,14 @@
+mtcars = dataset("datasets", "mtcars")
+
+formulae = [
+    @formula(MPG ~ HP + WT + Gear),
+    @formula(MPG ~ 0 + HP + WT + Gear),
+    @formula(MPG ~ HP + HP^2 + WT + WT * HP),
+    @formula(log(MPG) ~ log(HP) + log(WT))
+]
+
+for f in formulae
+    crrao_model = fit(f, mtcars, LinearRegression())
+    glm_model = lm(f, mtcars)
+    compare_models(crrao_model, glm_model, mtcars)
+end

test/numerical/frequentist/LogisticRegression.jl

@@ -0,0 +1,23 @@
+turnout = dataset("Zelig", "turnout")[1:100,:] # Take a subset of rows to reduce input size
+
+formulae = [
+    @formula(Vote ~ Age + Race + Income + Educate),
+    @formula(Vote ~ 0 + Age + Race + Income + Educate),
+    @formula(Vote ~ Age + Age^2 + Race + Income * Educate),
+    @formula(Vote ~ log(Age) + log(Educate))
+]
+
+links = [
+    (Logit(), LogitLink()),
+    (Probit(), ProbitLink()),
+    (Cloglog(), CloglogLink()),
+    (Cauchit(), CauchitLink())
+]
+
+for f in formulae
+    for (crrao_link, glm_link) in links
+        crrao_model = fit(f, turnout, LogisticRegression(), crrao_link)
+        glm_model = glm(f, turnout, Binomial(), glm_link)
+        compare_models(crrao_model, glm_model, turnout)
+    end
+end

test/numerical/frequentist/NegBinomialRegression.jl

@@ -0,0 +1,14 @@
+sanction = dataset("Zelig", "sanction")
+
+formulae = [
+    @formula(Num ~ Target + Coop + NCost),
+    @formula(Num ~ 0 + Target + Coop + NCost),
+    @formula(Num ~ Target + Target^2 + Coop + Coop * Target),
+    @formula(Num ~ log(Target) + log(Coop))
+]
+
+for f in formulae
+    crrao_model = fit(f, sanction, NegBinomRegression())
+    glm_model = glm(f, sanction, NegativeBinomial(), LogLink())
+    compare_models(crrao_model, glm_model, sanction)
+end

test/numerical/frequentist/PoissonRegression.jl

@@ -0,0 +1,14 @@
+sanction = dataset("Zelig", "sanction")
+
+formulae = [
+    @formula(Num ~ Target + Coop + NCost),
+    @formula(Num ~ 0 + Target + Coop + NCost),
+    @formula(Num ~ Target + Target^2 + Coop + Coop * Target),
+    @formula(Num ~ log(Target) + log(Coop))
+]
+
+for f in formulae
+    crrao_model = fit(f, sanction, PoissonRegression())
+    glm_model = glm(f, sanction, Poisson(), LogLink())
+    compare_models(crrao_model, glm_model, sanction)
+end

test/numerical/frequentist/tests.jl

@@ -0,0 +1,20 @@
+function compare_models(crrao_model, glm_model, df)
+    @test isapprox(coeftable(crrao_model).cols, coeftable(glm_model).cols)
+    @test isapprox(predict(crrao_model, df), predict(glm_model, df))
+end
+
+@testset "Linear Regression" begin
+    include("LinearRegression.jl")
+end
+
+@testset "Logistic Regression" begin
+    include("LogisticRegression.jl")
+end
+
+@testset "Poisson Regression" begin
+    include("PoissonRegression.jl")
+end
+
+@testset "Negative Binomial Regression" begin
+    include("NegBinomialRegression.jl")
+end

test/runtests.jl

@@ -1,4 +1,4 @@
-using CRRao, Test, StableRNGs, Logging, RDatasets, StatsModels
+using CRRao, Test, StableRNGs, Logging, RDatasets, StatsModels, GLM
 
 Logging.disable_logging(Logging.Warn)
 
@@ -23,4 +23,10 @@ CRRao.set_rng(StableRNG(123))
             include("basic/NegBinomialRegression.jl")
         end
     end
+
+    @testset "Numerical Tests" begin
+        @testset "Frequentist" begin
+            include("numerical/frequentist/tests.jl")
+        end
+    end
 end