contra contrasts

JuliaStats · Apr 12, 2023 · 386ede3 · 386ede3
1 parent 6be7c21
commit 386ede3
Show file tree

Hide file tree

Showing 4 changed files with 17 additions and 12 deletions.
diff --git a/test/bootstrap.jl b/test/bootstrap.jl
@@ -120,7 +120,8 @@ end
         contra = dataset(:contra)
         # need a model with fast=false to test that we only
         # copy the optimizer constraints for θ and not β
-        gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(), fast=false, progress=false)
+        gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli();
+                 fast=false, progress=false, contrasts=CONTRASTS[:contra])
         bs = parametricbootstrap(StableRNG(42), 100, gm0; hide_progress=true)
         # make sure we're not copying
         @test length(bs.lowerbd) == length(gm0.θ)

diff --git a/test/likelihoodratiotest.jl b/test/likelihoodratiotest.jl
@@ -92,8 +92,9 @@ end
     cc.usenum = ifelse.(cc.use .== "Y", 1 , 0)
     gmf = glm(@formula(usenum ~ 1+age+urban+livch), cc, Bernoulli());
     gmf2 = glm(@formula(usenum ~ 1+age+abs2(age)+urban+livch), cc, Bernoulli());
-    gm0 = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true, progress=false);
-    gm1 = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true, progress=false);
+    contrasts = CONTRASTS[:contra]
+    gm0 = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(); fast=true, progress=false, contrasts)
+    gm1 = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), contra, Bernoulli(); fast=true, progress=false, contrasts)
 
     lrt = likelihoodratiotest(gmf, gm1)
     @test [-2 * loglikelihood(gmf), deviance(gm1)] ≈ lrt.deviance
@@ -112,12 +113,12 @@ end
     @test length(lrt.formulae) == 2
 
     # mismatched links
-    gm_probit = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), ProbitLink(), fast=true, progress=false);
+    gm_probit = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), ProbitLink(); fast=true, progress=false, contrasts);
     @test_throws ArgumentError likelihoodratiotest(gmf, gm_probit)
     @test_throws ArgumentError likelihoodratiotest(gm0, gm_probit)
 
     # mismatched families
-    gm_poisson = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Poisson(), fast=true, progress=false);
+    gm_poisson = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Poisson(); fast=true, progress=false, contrasts);
     @test_throws ArgumentError likelihoodratiotest(gmf, gm_poisson)
     @test_throws ArgumentError likelihoodratiotest(gm0, gm_poisson)
 

diff --git a/test/pirls.jl b/test/pirls.jl
@@ -14,14 +14,16 @@ include("modelcache.jl")
 @testset "GLMM from MixedModel" begin
     f = first(gfms[:contra])
     d = dataset(:contra)
-    @test MixedModel(f, d, Bernoulli()) isa GeneralizedLinearMixedModel
-    @test MixedModel(f, d, Bernoulli(), ProbitLink()) isa GeneralizedLinearMixedModel
+    contrasts = CONTRASTS[:contra]
+    @test MixedModel(f, d, Bernoulli(); contrasts) isa GeneralizedLinearMixedModel
+    @test MixedModel(f, d, Bernoulli(), ProbitLink(); contrasts) isa GeneralizedLinearMixedModel
 end
 
 @testset "contra" begin
     contra = dataset(:contra)
+    contrasts = CONTRASTS[:contra]
     thin = 5
-    gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin)
+    gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin, contrasts)
     fitlog = gm0.optsum.fitlog
     @test length(fitlog) == (div(gm0.optsum.feval, thin) + 1) # for the initial value
     @test first(fitlog) == (gm0.optsum.initial, gm0.optsum.finitial)
@@ -48,7 +50,7 @@ end
     @test Distribution(gm0) == Distribution(gm0.resp)
     @test Link(gm0) == Link(gm0.resp)
 
-    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); progress=false);
+    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); progress=false, contrasts);
     @test isapprox(gm1.θ, [0.573054], atol=0.005)
     @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
     @test isapprox(deviance(gm1), 2361.54575, rtol=0.00001)
@@ -58,7 +60,7 @@ end
     @test nobs(gm0) == 1934
     refit!(gm0; fast=true, nAGQ=7, progress=false)
     @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
-    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); nAGQ=7, progress=false)
+    gm1 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); nAGQ=7, progress=false, contrasts)
     @test isapprox(deviance(gm1), 2360.8760, atol=0.001)
     @test gm1.β == gm1.beta
     @test gm1.θ == gm1.theta
@@ -82,7 +84,7 @@ end
     show(IOBuffer(), gm1)
     show(IOBuffer(), BlockDescription(gm0))
 
-    gm_slope = fit(MixedModel, gfms[:contra][2], contra, Bernoulli(); progress=false);
+    gm_slope = fit(MixedModel, gfms[:contra][2], contra, Bernoulli(); progress=false, contrasts);
     @test !issingular(gm_slope)
     @test issingular(gm_slope, zeros(5))
 

diff --git a/test/predict.jl b/test/predict.jl
@@ -111,8 +111,9 @@ end
 
     @testset "GLMM" begin
         contra = dataset(:contra)
+        contrasts = CONTRASTS[:contra]
         for fast in [true, false]
-            gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast, progress=false)
+            gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast, progress=false, contrasts)
 
             @test_throws ArgumentError predict(gm0, contra; type=:doh)