Support rank deficiency in MixedModels.jl extension (#66)

* ignore version specific manifest

* support rank deficiency in MixedModels extension

* version bump

* YAS

* test that shows we're aware of weirdness and it's consistent

.gitignore

@@ -22,3 +22,4 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+Manifest-v*.toml

Project.toml

@@ -1,7 +1,7 @@
 name = "Effects"
 uuid = "8f03c58b-bd97-4933-a826-f71b64d2cca2"
 authors = ["Beacon Biosignals, Inc."]
-version = "1.1.2"
+version = "1.2.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -42,6 +42,7 @@ Statistics = "1.6"
 StatsAPI = "1.6"
 StatsBase = "0.33, 0.34"
 StatsModels = "0.6.23, 0.7"
+Suppressor = "0.2"
 Tables = "1"
 Test = "1"
 TestSetExtensions = "3"
@@ -62,9 +63,10 @@ StandardizedPredictors = "5064a6a7-f8c2-40e2-8bdc-797ec6f1ae18"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 Vcov = "ec2bfdc2-55df-4fc9-b9ae-4958c2cf2486"
 
 [targets]
-test = ["Aqua", "DataFrames", "DataStructures", "GLM", "InlineStrings", "MixedModels", "MultipleTesting", "RDatasets", "StableRNGs", "StatsBase", "StatsModels", "StandardizedPredictors", "Test", "TestSetExtensions", "Vcov"]
+test = ["Aqua", "DataFrames", "DataStructures", "GLM", "InlineStrings", "MixedModels", "MultipleTesting", "RDatasets", "StableRNGs", "StatsBase", "StatsModels", "StandardizedPredictors", "Suppressor", "Test", "TestSetExtensions", "Vcov"]

ext/EffectsMixedModelsExt.jl

@@ -1,9 +1,32 @@
 module EffectsMixedModelsExt
 
+using DataFrames
 using Effects
 using MixedModels
+using StatsBase
+
+using Effects: typify, _difference_method!, _responsename
+using LinearAlgebra: diag
+using StatsModels: modelcols
 using GLM: Link
 
 Effects._model_link(m::GeneralizedLinearMixedModel, ::AutoInvLink) = Link(m)
 
+function Effects.effects!(reference_grid::DataFrame, model::MixedModel;
+                          eff_col=nothing, err_col=:err, typical=mean, invlink=identity,
+                          vcov=StatsBase.vcov)
+    # right now this is written for a RegressionModel and implicitly assumes
+    # the existence of an appropriate formula method
+    form = formula(model)
+    form_typical = typify(reference_grid, form, modelmatrix(model); typical=typical)
+    piv = view(model.feterm.piv, 1:(model.feterm.rank))
+    X = view(modelcols(form_typical, reference_grid), :, piv)
+    eff = X * fixef(model)
+    err = sqrt.(diag(X * vcov(model)[piv, piv] * X'))
+    _difference_method!(eff, err, model, invlink)
+    reference_grid[!, something(eff_col, _responsename(model))] = eff
+    reference_grid[!, err_col] = err
+    return reference_grid
+end
+
 end # module

test/mixedmodels.jl

@@ -0,0 +1,35 @@
+using DataFrames
+using Effects
+using MixedModels
+using StableRNGs
+using Suppressor
+using Test
+
+rng = StableRNG(42)
+x = rand(rng, 100)
+data = (x=x, x2=1.5 .* x, y=rand(rng, [0, 1], 100), z=repeat('A':'T', 5))
+model = @suppress fit(MixedModel, @formula(y ~ x + x2 + (1 | z)), data; progress=false)
+dropped_idx = model.feterm.piv[end]
+dropped_coef = coefnames(model)[dropped_idx]
+kept_coef = last(fixefnames(model))
+
+# due to the vagaries of numerical linear algebra and pivoting, it's
+# not completely deterministic which of x and x2 gets dropped, though it
+# will tend to be x2
+design = Dict(Symbol(kept_coef) => [0])
+eff = effects(design, model)
+@test dropped_coef ∉ names(eff)
+@test kept_coef ∈ names(eff)
+@test !isnan(only(eff.err))
+@test eff.y - eff.err ≈ eff.lower
+@test eff.y + eff.err ≈ eff.upper
+
+# there is one bit of weirdness -- the removed coefficients behave like any other
+# variable in the "design" that is missing from the model.
+design = Dict(Symbol(dropped_coef) => [0])
+eff_dropped = effects(design, model)
+
+design = Dict(:q => [0])
+eff_not_present = effects(design, model)
+
+@test Matrix(eff_dropped) ≈ Matrix(eff_not_present)

test/runtests.jl

@@ -23,4 +23,10 @@ using TestSetExtensions
     @testset "emmeans" begin
         include("emmeans.jl")
     end
+
+    @static if VERSION >= v"1.9"
+        @testset "MixedModels.jl" begin
+            include("mixedmodels.jl")
+        end
+    end
 end