From e1eb68965c30baf14a47edb72b94fdd297e557ba Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 16:55:03 -0500 Subject: [PATCH 01/15] GLM ext --- Project.toml | 12 +++++- ext/EffectsGLMExt.jl | 25 +++++++++++ ext/EffectsMixedModelsExt.jl | 10 +++++ src/Effects.jl | 2 +- src/regressionmodel.jl | 52 ++++++++++++++++++++-- test/delta_method.jl | 83 +++++++++++++++++++----------------- 6 files changed, 141 insertions(+), 43 deletions(-) create mode 100644 ext/EffectsGLMExt.jl create mode 100644 ext/EffectsMixedModelsExt.jl diff --git a/Project.toml b/Project.toml index f31d38b..dda2b95 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Effects" uuid = "8f03c58b-bd97-4933-a826-f71b64d2cca2" authors = ["Beacon Biosignals, Inc."] -version = "1.0.2" +version = "1.1.0" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -10,10 +10,19 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +[weakdeps] +GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" +MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" + +[extensions] +EffectsGLMExt = "GLM" +EffectsMixedModelsExt = ["GLM", "MixedModels"] + [compat] Aqua = "0.5, 0.6" Combinatorics = "1" @@ -27,6 +36,7 @@ MultipleTesting = "0.5, 0.6" RDatasets = "0.7.7" StableRNGs = "1.0" StandardizedPredictors = "0.1, 1" +StatsAPI = "1.6" StatsBase = "0.33, 0.34" StatsModels = "0.6.23, 0.7" Tables = "1" diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl new file mode 100644 index 0000000..e19cb51 --- /dev/null +++ b/ext/EffectsGLMExt.jl @@ -0,0 +1,25 @@ +module EffectsGLMExt + +using Effects + +using GLM: AbstractGLM, Link, mueta, linkinv +using StatsAPI: RegressionModel +using StatsModels: TableRegressionModel + +# we keep RegressionModel so that this will also be used +# for MixedModels.jl + +Effects._link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) + +function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, + model::RegressionModel, + ::AutoInvLink) where {T <: AbstractFloat} + link = Effects._link(model) + err .*= mueta.(link, eff) + eff .= linkinv.(link, eff) + + return err +end + + +end # module diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl new file mode 100644 index 0000000..c125e00 --- /dev/null +++ b/ext/EffectsMixedModelsExt.jl @@ -0,0 +1,10 @@ +module EffectsMixedModelsExt + +using Effects +using MixedModels + +using GLM: Link + +Effects._link(m::GeneralizedLinearMixedModel) = Link(m.glm) + +end # module diff --git a/src/Effects.jl b/src/Effects.jl index 76dafea..b9a1db5 100644 --- a/src/Effects.jl +++ b/src/Effects.jl @@ -15,7 +15,7 @@ using ForwardDiff include("typical.jl") include("regressionmodel.jl") -export effects, effects!, expand_grid +export effects, effects!, expand_grid, AutoInvLink include("emmeans.jl") export emmeans, empairs diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index 0b4da18..1622fd0 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -1,3 +1,25 @@ +# determine AutoInvLink automatically if a package with an appropriate extension +# is loaded + +""" + AutoInvLink + +Singleton type indicating that the inverse link should be automatically +determined from the model type. + +This will only work for model types with an appropriate extension on +Julia 1.9+. If an appropriate extension is not defined, then an error +will occur. + +Currently, this is only implemented for GLM.jl and MixedModels.jl +""" +struct AutoInvLink end + +# helper function to override in extensions +function _link(::Any) + throw(ArgumentError("No appropriate extension is loaded for automatic " * + "determination of the inverse link for this model type")) +end """ effects!(reference_grid::DataFrame, model::RegressionModel; eff_col=nothing, err_col=:err, typical=mean, invlink=identity, @@ -54,6 +76,10 @@ Pointwise standard errors are written into the column specified by `err_col`. automatic differentiation. This means that the `invlink` function must be differentiable and should not involve inplace operations. +The special singleton value `AutoInvLink()` can be used to specify that +the appropriate inverse link should be determined and, where possible, a +direct or analytic computation of the derivative is used. + Effects are computed using the model's variance-covariance matrix, which is computed by default using `StatsBas.vcov`. Alternative methods such as the sandwich estimator or robust estimators can be used by specifying `vcov`, @@ -67,7 +93,6 @@ using Vcov myvcov(x) = Base.Fix2(vcov, Vcov.robust()) ``` - The reference grid must contain columns for all predictors in the formula. (Interactions are computed automatically.) Contrasts must match the contrasts used to fit the model; using other contrasts will lead to undefined behavior. @@ -87,6 +112,8 @@ The approach for computing effect is based on the effects plots described here: Fox, John (2003). Effect Displays in R for Generalised Linear Models. Journal of Statistical Software. Vol. 8, No. 15 + +See also [`AutoInvLink`](@ref). """ function effects!(reference_grid::DataFrame, model::RegressionModel; eff_col=nothing, err_col=:err, typical=mean, invlink=identity, @@ -99,8 +126,7 @@ function effects!(reference_grid::DataFrame, model::RegressionModel; eff = X * coef(model) err = sqrt.(diag(X * vcov(model) * X')) if invlink !== identity - err .*= ForwardDiff.derivative.(invlink, eff) - eff .= invlink.(eff) + _difference_method!(eff, err, model, invlink) end reference_grid[!, something(eff_col, _responsename(model))] = eff reference_grid[!, err_col] = err @@ -110,6 +136,24 @@ function effects!(reference_grid::DataFrame, model::RegressionModel; # return (; reference_grid..., depvar => eff, err_col => err) end +# TODO: support the transformation method +# in addition to the difference method +# xref https://github.com/JuliaStats/GLM.jl/blob/c13577eaf3f418c58020534dd407532ee57f219b/src/glmfit.jl#L773-L783 + +function _difference_method!(eff::Vector{T}, err::Vector{T}, + ::RegressionModel, + invlink) where {T <: AbstractFloat} + + err .*= ForwardDiff.derivative.(invlink, eff) + eff .= invlink.(eff) + return eff, err +end + +function ForwardDiff.derivative(::AutoInvLink, ::Real) + throw(ArgumentError("No appropriate extension is loaded for automatic " * + "determination of the inverse link for this model type")) +end + """ expand_grid(design) @@ -155,6 +199,8 @@ fully-crossed design. Additionally, two extra columns are created representing the lower and upper edges of the confidence interval. The default `level=nothing` corresponds to [resp-err, resp+err], while `level=0.95` corresponds to the 95% confidence interval. + +See also [`AutoInvLink`](@ref). """ function effects(design::AbstractDict, model::RegressionModel; eff_col=nothing, err_col=:err, typical=mean, diff --git a/test/delta_method.jl b/test/delta_method.jl index b8960d4..c85f953 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -8,9 +8,9 @@ using Test @testset "transformed response" begin dat = rdataset("car", "Prestige") - model = lm(@formula(log(Prestige) ~ 1 + Income * Education), dat) design = Dict(:Income => [1, 2], :Education => [3, 4]) + model = lm(@formula(log(Prestige) ~ 1 + Income * Education), dat) eff_original_scale = effects(design, model; invlink=exp) eff_logscale = effects(design, model) @test all(eff_logscale[!, 3] .≈ log.(eff_original_scale[!, 3])) @@ -37,55 +37,62 @@ using Test @test isapprox(only(eff_emm_trans.err), 1.07; atol=0.005) @test isapprox(only(eff_emm_trans.lower), 45.3; atol=0.05) @test isapprox(only(eff_emm_trans.upper), 47.5; atol=0.05) + + @testset "AutoInvLink doesn't work for models without a Link()" begin + # this goes through _link(::Any) + @test_throws ArgumentError effects(design, model; invlink=AutoInvLink()) + end end @testset "link function" begin dat = rdataset("car", "Cowles") dat[!, :vol] = dat.Volunteer .== "yes" model = glm(@formula(vol ~ Extraversion * Neuroticism), dat, Bernoulli()) - invlink = Base.Fix1(GLM.linkinv, Link(model.model)) design = Dict(:Extraversion => [13], - :Neuroticism => [16]) + :Neuroticism => [16]) X = [1.0 13.0 16.0 13 * 16] + iv = Base.Fix1(GLM.linkinv, Link(model.model)) + invlinks = [iv, AutoInvLink()] + @testset "invlink = $invlink" for invlink in invlinks + for level in [0.68, 0.95] + eff = effects(design, model; invlink, level) - for level in [0.68, 0.95] - eff = effects(design, model; invlink, level) + # compare with results from GLM.predict + pred = DataFrame(predict(model.model, X; + interval=:confidence, + interval_method=:delta, + level)) + @test all(pred.prediction .≈ eff.vol) + @test all(isapprox.(pred.lower, eff.lower; atol=0.001)) + @test all(isapprox.(pred.upper, eff.upper; atol=0.001)) - # compare with results from GLM.predict - pred = DataFrame(predict(model.model, X; - interval=:confidence, - interval_method=:delta, - level)) - @test all(pred.prediction .≈ eff.vol) - @test all(isapprox.(pred.lower, eff.lower; atol=0.001)) - @test all(isapprox.(pred.upper, eff.upper; atol=0.001)) + eff_trans = effects(design, model; level) + transform!(eff_trans, + :vol => ByRow(iv), + :lower => ByRow(iv), + :upper => ByRow(iv); renamecols=false) + # for this model, things play out nicely + @test all(eff_trans.vol .≈ eff.vol) + @test all(isapprox.(eff_trans.lower, eff.lower; atol=0.001)) + @test all(isapprox.(eff_trans.upper, eff.upper; atol=0.001)) + end - eff_trans = effects(design, model; level) - transform!(eff_trans, - :vol => ByRow(invlink), - :lower => ByRow(invlink), - :upper => ByRow(invlink); renamecols=false) - # for this model, things play out nicely - @test all(eff_trans.vol .≈ eff.vol) - @test all(isapprox.(eff_trans.lower, eff.lower; atol=0.001)) - @test all(isapprox.(eff_trans.upper, eff.upper; atol=0.001)) - end + # compare with results from emmeans in R + # emmeans(model, ~ neuroticism * extraversion, level=0.68) + eff_emm = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model) + @test isapprox(only(eff_emm.vol), -0.347; atol=0.005) + @test isapprox(only(eff_emm.err), 0.0549; atol=0.005) + @test isapprox(only(eff_emm.lower), -0.402; atol=0.005) + @test isapprox(only(eff_emm.upper), -0.292; atol=0.005) - # compare with results from emmeans in R - # emmeans(model, ~ neuroticism * extraversion, level=0.68) - eff_emm = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model) - @test isapprox(only(eff_emm.vol), -0.347; atol=0.005) - @test isapprox(only(eff_emm.err), 0.0549; atol=0.005) - @test isapprox(only(eff_emm.lower), -0.402; atol=0.005) - @test isapprox(only(eff_emm.upper), -0.292; atol=0.005) - - # emmeans(model, ~ neuroticism * extraversion, level=0.68, transform="response") - eff_emm_trans = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model; - invlink) - @test isapprox(only(eff_emm_trans.vol), 0.414; atol=0.005) - @test isapprox(only(eff_emm_trans.err), 0.0133; atol=0.005) - @test isapprox(only(eff_emm_trans.lower), 0.401; atol=0.005) - @test isapprox(only(eff_emm_trans.upper), 0.427; atol=0.005) + # emmeans(model, ~ neuroticism * extraversion, level=0.68, transform="response") + eff_emm_trans = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model; + invlink) + @test isapprox(only(eff_emm_trans.vol), 0.414; atol=0.005) + @test isapprox(only(eff_emm_trans.err), 0.0133; atol=0.005) + @test isapprox(only(eff_emm_trans.lower), 0.401; atol=0.005) + @test isapprox(only(eff_emm_trans.upper), 0.427; atol=0.005) + end end @testset "identity by another name" begin From 6726f88edc47ad277aa76c854b707b30b74b060e Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 16:57:16 -0500 Subject: [PATCH 02/15] better --- ext/EffectsGLMExt.jl | 7 ++++--- ext/EffectsMixedModelsExt.jl | 11 ++++++++--- src/regressionmodel.jl | 9 +++------ test/delta_method.jl | 3 +-- 4 files changed, 16 insertions(+), 14 deletions(-) diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index e19cb51..e76bcfb 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -9,12 +9,13 @@ using StatsModels: TableRegressionModel # we keep RegressionModel so that this will also be used # for MixedModels.jl -Effects._link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) +_link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) +_link(m::AbstractGLM) = Link(m) function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, - model::RegressionModel, + model::Union{TableRegressionModel{<:AbstractGLM}, AbstractGLM}, ::AutoInvLink) where {T <: AbstractFloat} - link = Effects._link(model) + link = _link(model) err .*= mueta.(link, eff) eff .= linkinv.(link, eff) diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl index c125e00..fa0faec 100644 --- a/ext/EffectsMixedModelsExt.jl +++ b/ext/EffectsMixedModelsExt.jl @@ -3,8 +3,13 @@ module EffectsMixedModelsExt using Effects using MixedModels -using GLM: Link - -Effects._link(m::GeneralizedLinearMixedModel) = Link(m.glm) +function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, + model::GeneralizedLinearMixedModel, + ::AutoInvLink) where {T <: AbstractFloat} + link = Link(model) + err .*= mueta.(link, eff) + eff .= linkinv.(link, eff) + return err +end end # module diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index 1622fd0..31ab9f0 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -15,11 +15,6 @@ Currently, this is only implemented for GLM.jl and MixedModels.jl """ struct AutoInvLink end -# helper function to override in extensions -function _link(::Any) - throw(ArgumentError("No appropriate extension is loaded for automatic " * - "determination of the inverse link for this model type")) -end """ effects!(reference_grid::DataFrame, model::RegressionModel; eff_col=nothing, err_col=:err, typical=mean, invlink=identity, @@ -149,7 +144,9 @@ function _difference_method!(eff::Vector{T}, err::Vector{T}, return eff, err end -function ForwardDiff.derivative(::AutoInvLink, ::Real) +function _difference_method!(::Vector{T}, ::Vector{T}, + ::RegressionModel, + ::AutoInvLink) where {T <: AbstractFloat} throw(ArgumentError("No appropriate extension is loaded for automatic " * "determination of the inverse link for this model type")) end diff --git a/test/delta_method.jl b/test/delta_method.jl index c85f953..1734a38 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -38,8 +38,7 @@ using Test @test isapprox(only(eff_emm_trans.lower), 45.3; atol=0.05) @test isapprox(only(eff_emm_trans.upper), 47.5; atol=0.05) - @testset "AutoInvLink doesn't work for models without a Link()" begin - # this goes through _link(::Any) + @testset "AutoInvLink fails gracefully" begin @test_throws ArgumentError effects(design, model; invlink=AutoInvLink()) end end From ea08c3d1a963574b6350964dc174ddb34ce4cbac Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 17:11:07 -0500 Subject: [PATCH 03/15] test mixed models extension --- ext/EffectsMixedModelsExt.jl | 1 + test/delta_method.jl | 14 ++++++++++++++ 2 files changed, 15 insertions(+) diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl index fa0faec..69e978a 100644 --- a/ext/EffectsMixedModelsExt.jl +++ b/ext/EffectsMixedModelsExt.jl @@ -2,6 +2,7 @@ module EffectsMixedModelsExt using Effects using MixedModels +using GLM: Link, mueta, linkinv function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, model::GeneralizedLinearMixedModel, diff --git a/test/delta_method.jl b/test/delta_method.jl index 1734a38..417d419 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -2,6 +2,7 @@ using DataFrames using Effects using GLM using LinearAlgebra +using MixedModels using RDatasets: dataset as rdataset using StableRNGs using Test @@ -94,6 +95,19 @@ end end end +@testset "link function in a MixedModel" begin + model = fit(MixedModel, + @formula(use ~ 1 + age + (1|urban)), + MixedModels.dataset(:contra), + Bernoulli(); progress=false) + design = Dict(:age => -10:10) + eff_manual = effects(design, model; + invlink=Base.Fix1(GLM.linkinv, Link(model))) + eff_auto = effects(design, model; invlink=AutoInvLink()) + + @test all(isapprox.(Matrix(eff_manual), Matrix(eff_auto))) +end + @testset "identity by another name" begin b0, b1, b2, b1_2 = beta = [0.0, 1.0, 1.0, -1.0] x = collect(-10:10) From 9a447a9718032d7956941397dedcce3ada722cc5 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 17:24:32 -0500 Subject: [PATCH 04/15] better test support for Julia <1.9 --- src/regressionmodel.jl | 3 +++ test/delta_method.jl | 30 +++++++++++++++++++----------- test/runtests.jl | 2 +- 3 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index 31ab9f0..e1e8a57 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -147,6 +147,9 @@ end function _difference_method!(::Vector{T}, ::Vector{T}, ::RegressionModel, ::AutoInvLink) where {T <: AbstractFloat} + @static if VERSION < v"1.9" + @error "AutoInvLink requires extensions and is thus not available on Julia < 1.9." + end throw(ArgumentError("No appropriate extension is loaded for automatic " * "determination of the inverse link for this model type")) end diff --git a/test/delta_method.jl b/test/delta_method.jl index 417d419..6d61710 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -40,6 +40,8 @@ using Test @test isapprox(only(eff_emm_trans.upper), 47.5; atol=0.05) @testset "AutoInvLink fails gracefully" begin + # this should work even pre Julia 1.9 because by definition + # no extension is loaded @test_throws ArgumentError effects(design, model; invlink=AutoInvLink()) end end @@ -52,7 +54,11 @@ end :Neuroticism => [16]) X = [1.0 13.0 16.0 13 * 16] iv = Base.Fix1(GLM.linkinv, Link(model.model)) - invlinks = [iv, AutoInvLink()] + @static if VERSION >= v"1.9" + invlinks = [iv, AutoInvLink()] + else + invlinks = [iv] + end @testset "invlink = $invlink" for invlink in invlinks for level in [0.68, 0.95] eff = effects(design, model; invlink, level) @@ -95,17 +101,19 @@ end end end -@testset "link function in a MixedModel" begin - model = fit(MixedModel, - @formula(use ~ 1 + age + (1|urban)), - MixedModels.dataset(:contra), - Bernoulli(); progress=false) - design = Dict(:age => -10:10) - eff_manual = effects(design, model; - invlink=Base.Fix1(GLM.linkinv, Link(model))) - eff_auto = effects(design, model; invlink=AutoInvLink()) +@static if VERSION >= v"1.9" + @testset "link function in a MixedModel" begin + model = fit(MixedModel, + @formula(use ~ 1 + age + (1|urban)), + MixedModels.dataset(:contra), + Bernoulli(); progress=false) + design = Dict(:age => -10:10) + eff_manual = effects(design, model; + invlink=Base.Fix1(GLM.linkinv, Link(model))) + eff_auto = effects(design, model; invlink=AutoInvLink()) - @test all(isapprox.(Matrix(eff_manual), Matrix(eff_auto))) + @test all(isapprox.(Matrix(eff_manual), Matrix(eff_auto))) + end end @testset "identity by another name" begin diff --git a/test/runtests.jl b/test/runtests.jl index 6dbc1de..94d6117 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using Effects using Test @testset "Aqua" begin - Aqua.test_all(Effects; ambiguities=false) + Aqua.test_all(Effects; ambiguities=false, project_toml_formatting) end @testset "TypicalTerm" begin From 1433b1e1ffa7614bc55b1bad5c1a28b035375fb6 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 17:27:12 -0500 Subject: [PATCH 05/15] YAS --- .github/workflows/YASG.yml | 1 + ext/EffectsGLMExt.jl | 8 ++++---- ext/EffectsMixedModelsExt.jl | 6 +++--- format/run.jl | 2 +- src/regressionmodel.jl | 15 +++++++-------- test/delta_method.jl | 27 ++++++++++++++------------- 6 files changed, 30 insertions(+), 29 deletions(-) diff --git a/.github/workflows/YASG.yml b/.github/workflows/YASG.yml index 4896c32..5aee7b9 100644 --- a/.github/workflows/YASG.yml +++ b/.github/workflows/YASG.yml @@ -11,6 +11,7 @@ on: - 'src/**' - 'test/**' - 'docs/**' + - 'ext/**' - '.github/workflows/YASG.yml' - 'format/**' jobs: diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index e76bcfb..52bc63b 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -12,9 +12,10 @@ using StatsModels: TableRegressionModel _link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) _link(m::AbstractGLM) = Link(m) -function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, - model::Union{TableRegressionModel{<:AbstractGLM}, AbstractGLM}, - ::AutoInvLink) where {T <: AbstractFloat} +function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, + model::Union{TableRegressionModel{<:AbstractGLM}, + AbstractGLM}, + ::AutoInvLink) where {T<:AbstractFloat} link = _link(model) err .*= mueta.(link, eff) eff .= linkinv.(link, eff) @@ -22,5 +23,4 @@ function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, return err end - end # module diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl index 69e978a..7d3fe2a 100644 --- a/ext/EffectsMixedModelsExt.jl +++ b/ext/EffectsMixedModelsExt.jl @@ -4,9 +4,9 @@ using Effects using MixedModels using GLM: Link, mueta, linkinv -function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, - model::GeneralizedLinearMixedModel, - ::AutoInvLink) where {T <: AbstractFloat} +function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, + model::GeneralizedLinearMixedModel, + ::AutoInvLink) where {T<:AbstractFloat} link = Link(model) err .*= mueta.(link, eff) eff .= linkinv.(link, eff) diff --git a/format/run.jl b/format/run.jl index 5e7150d..87bfc57 100644 --- a/format/run.jl +++ b/format/run.jl @@ -3,7 +3,7 @@ using JuliaFormatter function main() perfect = true # note: keep in sync with `.github/workflows/format-check.yml` - for d in ["src/", "test/", "docs/"] + for d in ["src/", "ext/", "test/", "docs/"] @info "...linting $d ..." dir_perfect = format(d, YASStyle()) perfect = perfect && dir_perfect diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index e1e8a57..e63ae19 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -135,22 +135,21 @@ end # in addition to the difference method # xref https://github.com/JuliaStats/GLM.jl/blob/c13577eaf3f418c58020534dd407532ee57f219b/src/glmfit.jl#L773-L783 -function _difference_method!(eff::Vector{T}, err::Vector{T}, - ::RegressionModel, - invlink) where {T <: AbstractFloat} - +function _difference_method!(eff::Vector{T}, err::Vector{T}, + ::RegressionModel, + invlink) where {T<:AbstractFloat} err .*= ForwardDiff.derivative.(invlink, eff) eff .= invlink.(eff) return eff, err end -function _difference_method!(::Vector{T}, ::Vector{T}, - ::RegressionModel, - ::AutoInvLink) where {T <: AbstractFloat} +function _difference_method!(::Vector{T}, ::Vector{T}, + ::RegressionModel, + ::AutoInvLink) where {T<:AbstractFloat} @static if VERSION < v"1.9" @error "AutoInvLink requires extensions and is thus not available on Julia < 1.9." end - throw(ArgumentError("No appropriate extension is loaded for automatic " * + throw(ArgumentError("No appropriate extension is loaded for automatic " * "determination of the inverse link for this model type")) end diff --git a/test/delta_method.jl b/test/delta_method.jl index 6d61710..7fa3c19 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -51,7 +51,7 @@ end dat[!, :vol] = dat.Volunteer .== "yes" model = glm(@formula(vol ~ Extraversion * Neuroticism), dat, Bernoulli()) design = Dict(:Extraversion => [13], - :Neuroticism => [16]) + :Neuroticism => [16]) X = [1.0 13.0 16.0 13 * 16] iv = Base.Fix1(GLM.linkinv, Link(model.model)) @static if VERSION >= v"1.9" @@ -65,18 +65,18 @@ end # compare with results from GLM.predict pred = DataFrame(predict(model.model, X; - interval=:confidence, - interval_method=:delta, - level)) + interval=:confidence, + interval_method=:delta, + level)) @test all(pred.prediction .≈ eff.vol) @test all(isapprox.(pred.lower, eff.lower; atol=0.001)) @test all(isapprox.(pred.upper, eff.upper; atol=0.001)) eff_trans = effects(design, model; level) transform!(eff_trans, - :vol => ByRow(iv), - :lower => ByRow(iv), - :upper => ByRow(iv); renamecols=false) + :vol => ByRow(iv), + :lower => ByRow(iv), + :upper => ByRow(iv); renamecols=false) # for this model, things play out nicely @test all(eff_trans.vol .≈ eff.vol) @test all(isapprox.(eff_trans.lower, eff.lower; atol=0.001)) @@ -92,7 +92,8 @@ end @test isapprox(only(eff_emm.upper), -0.292; atol=0.005) # emmeans(model, ~ neuroticism * extraversion, level=0.68, transform="response") - eff_emm_trans = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model; + eff_emm_trans = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), + model; invlink) @test isapprox(only(eff_emm_trans.vol), 0.414; atol=0.005) @test isapprox(only(eff_emm_trans.err), 0.0133; atol=0.005) @@ -103,13 +104,13 @@ end @static if VERSION >= v"1.9" @testset "link function in a MixedModel" begin - model = fit(MixedModel, - @formula(use ~ 1 + age + (1|urban)), - MixedModels.dataset(:contra), + model = fit(MixedModel, + @formula(use ~ 1 + age + (1 | urban)), + MixedModels.dataset(:contra), Bernoulli(); progress=false) design = Dict(:age => -10:10) - eff_manual = effects(design, model; - invlink=Base.Fix1(GLM.linkinv, Link(model))) + eff_manual = effects(design, model; + invlink=Base.Fix1(GLM.linkinv, Link(model))) eff_auto = effects(design, model; invlink=AutoInvLink()) @test all(isapprox.(Matrix(eff_manual), Matrix(eff_auto))) From b0355e49a7461ff31234c6c8e7e837a42734a433 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 17:30:50 -0500 Subject: [PATCH 06/15] remove cruft --- ext/EffectsGLMExt.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index 52bc63b..dd9a746 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -6,9 +6,6 @@ using GLM: AbstractGLM, Link, mueta, linkinv using StatsAPI: RegressionModel using StatsModels: TableRegressionModel -# we keep RegressionModel so that this will also be used -# for MixedModels.jl - _link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) _link(m::AbstractGLM) = Link(m) From 6ac86f846d205b9dd5ab6dbb09d18be1baae6fc9 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 17:37:32 -0500 Subject: [PATCH 07/15] oops --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 94d6117..6dbc1de 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using Effects using Test @testset "Aqua" begin - Aqua.test_all(Effects; ambiguities=false, project_toml_formatting) + Aqua.test_all(Effects; ambiguities=false) end @testset "TypicalTerm" begin From 178f93ea5e7957fb64af06d7fd2669ebe544f469 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 18:12:47 -0500 Subject: [PATCH 08/15] disable Project.toml formatting check on LTS --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6dbc1de..c40fb4a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,9 @@ using Effects using Test @testset "Aqua" begin - Aqua.test_all(Effects; ambiguities=false) + # layout of weakdeps and extensions etc. differs between pre 1.9 and 1.9+ + project_toml_formatting = VERSION >= v"1.9" + Aqua.test_all(Effects; ambiguities=false, project_toml_formatting) end @testset "TypicalTerm" begin From 66459e7838075b1605e9f696133dd5305a5aee8c Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 19 Aug 2023 20:16:56 -0500 Subject: [PATCH 09/15] remove unused method --- ext/EffectsGLMExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index dd9a746..0fdd664 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -6,8 +6,8 @@ using GLM: AbstractGLM, Link, mueta, linkinv using StatsAPI: RegressionModel using StatsModels: TableRegressionModel +# TODO: upstream a Link(::TableRegressionModel{<:AbstractGLM}) _link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) -_link(m::AbstractGLM) = Link(m) function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, model::Union{TableRegressionModel{<:AbstractGLM}, From b6cf8e2d4079e91100a0761ef43c55b239641bfd Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 13:35:17 -0500 Subject: [PATCH 10/15] reduce duplication Co-authored-by: Alex Arslan --- ext/EffectsGLMExt.jl | 19 +++++-------------- ext/EffectsMixedModelsExt.jl | 11 ++--------- src/regressionmodel.jl | 34 ++++++++++++++++++++-------------- test/delta_method.jl | 1 + 4 files changed, 28 insertions(+), 37 deletions(-) diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index 0fdd664..588aa02 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -2,22 +2,13 @@ module EffectsGLMExt using Effects -using GLM: AbstractGLM, Link, mueta, linkinv -using StatsAPI: RegressionModel +using GLM: AbstractGLM, Link, Link01, inverselink using StatsModels: TableRegressionModel # TODO: upstream a Link(::TableRegressionModel{<:AbstractGLM}) -_link(m::TableRegressionModel{<:AbstractGLM}) = Link(m.model) - -function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, - model::Union{TableRegressionModel{<:AbstractGLM}, - AbstractGLM}, - ::AutoInvLink) where {T<:AbstractFloat} - link = _link(model) - err .*= mueta.(link, eff) - eff .= linkinv.(link, eff) - - return err -end +Effects._model_link(m::TableRegressionModel{<:AbstractGLM}, ::AutoInvLink) = Link(m.model) +Effects._model_link(m::AbstractGLM, ::AutoInvLink) = Link(m) +Effects._invlink_and_deriv(link::Link01, η) = inverselink(link, η)[1:2:3] # (µ, 1 - µ, dμdη) +Effects._invlink_and_deriv(link::Link, η) = inverselink(link, η)[1:2] # (µ, dμdη, NaN) end # module diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl index 7d3fe2a..9d77716 100644 --- a/ext/EffectsMixedModelsExt.jl +++ b/ext/EffectsMixedModelsExt.jl @@ -2,15 +2,8 @@ module EffectsMixedModelsExt using Effects using MixedModels -using GLM: Link, mueta, linkinv +using GLM: Link -function Effects._difference_method!(eff::Vector{T}, err::Vector{T}, - model::GeneralizedLinearMixedModel, - ::AutoInvLink) where {T<:AbstractFloat} - link = Link(model) - err .*= mueta.(link, eff) - eff .= linkinv.(link, eff) - return err -end +Effects._model_link(m::GeneralizedLinearMixedModel, ::AutoInvLink) = Link(m) end # module diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index e63ae19..cb72bbb 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -120,9 +120,7 @@ function effects!(reference_grid::DataFrame, model::RegressionModel; X = modelcols(form_typical, reference_grid) eff = X * coef(model) err = sqrt.(diag(X * vcov(model) * X')) - if invlink !== identity - _difference_method!(eff, err, model, invlink) - end + _difference_method!(eff, err, model, invlink) reference_grid[!, something(eff_col, _responsename(model))] = eff reference_grid[!, err_col] = err return reference_grid @@ -135,17 +133,13 @@ end # in addition to the difference method # xref https://github.com/JuliaStats/GLM.jl/blob/c13577eaf3f418c58020534dd407532ee57f219b/src/glmfit.jl#L773-L783 -function _difference_method!(eff::Vector{T}, err::Vector{T}, - ::RegressionModel, - invlink) where {T<:AbstractFloat} - err .*= ForwardDiff.derivative.(invlink, eff) - eff .= invlink.(eff) - return eff, err -end - -function _difference_method!(::Vector{T}, ::Vector{T}, - ::RegressionModel, - ::AutoInvLink) where {T<:AbstractFloat} +_invlink_and_deriv(invlink, η) = (invlink(η), ForwardDiff.derivative(invlink, η)) +_invlink_and_deriv(::typeof(identity), η) = (η, 1) +# this isn't the best name because it sometimes returns the inverse link and sometimes the link (Link()) +# for now, this is private API, but we should see how this goes and whether we can make it public API +# so local extensions (instead of Package-Extensions) are better supported +_model_link(::RegressionModel, invlink::Function) = invlink +function _model_link(::RegressionModel, ::AutoInvLink) @static if VERSION < v"1.9" @error "AutoInvLink requires extensions and is thus not available on Julia < 1.9." end @@ -153,6 +147,18 @@ function _difference_method!(::Vector{T}, ::Vector{T}, "determination of the inverse link for this model type")) end +function _difference_method!(eff::Vector{T}, err::Vector{T}, + m::RegressionModel, + invlink) where {T<:AbstractFloat} + link = _model_link(m, invlink) + @inbounds for i in eachindex(eff, err) + μ, dμdη = _invlink_and_deriv(link, eff[i]) + err[i] *= dμdη + eff[i] = μ + end + return eff, err +end + """ expand_grid(design) diff --git a/test/delta_method.jl b/test/delta_method.jl index 7fa3c19..5185252 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -56,6 +56,7 @@ end iv = Base.Fix1(GLM.linkinv, Link(model.model)) @static if VERSION >= v"1.9" invlinks = [iv, AutoInvLink()] + @test Effects._model_link(model, AutoInvLink()) == Effects._model_link(model.model, AutoInvLink()) else invlinks = [iv] end From e79867b7bc2c0fb941fcf23b9226a296e44dd343 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 18:38:32 +0000 Subject: [PATCH 11/15] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- ext/EffectsMixedModelsExt.jl | 2 +- src/regressionmodel.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/EffectsMixedModelsExt.jl b/ext/EffectsMixedModelsExt.jl index 9d77716..ff83274 100644 --- a/ext/EffectsMixedModelsExt.jl +++ b/ext/EffectsMixedModelsExt.jl @@ -4,6 +4,6 @@ using Effects using MixedModels using GLM: Link -Effects._model_link(m::GeneralizedLinearMixedModel, ::AutoInvLink) = Link(m) +Effects._model_link(m::GeneralizedLinearMixedModel, ::AutoInvLink) = Link(m) end # module diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index cb72bbb..0c02298 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -139,7 +139,7 @@ _invlink_and_deriv(::typeof(identity), η) = (η, 1) # for now, this is private API, but we should see how this goes and whether we can make it public API # so local extensions (instead of Package-Extensions) are better supported _model_link(::RegressionModel, invlink::Function) = invlink -function _model_link(::RegressionModel, ::AutoInvLink) +function _model_link(::RegressionModel, ::AutoInvLink) @static if VERSION < v"1.9" @error "AutoInvLink requires extensions and is thus not available on Julia < 1.9." end From b925065b955998898a82077b5b9aeb05ba2e634a Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 14:49:07 -0500 Subject: [PATCH 12/15] Update test/delta_method.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- test/delta_method.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/delta_method.jl b/test/delta_method.jl index 5185252..3a3e3f7 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -56,7 +56,8 @@ end iv = Base.Fix1(GLM.linkinv, Link(model.model)) @static if VERSION >= v"1.9" invlinks = [iv, AutoInvLink()] - @test Effects._model_link(model, AutoInvLink()) == Effects._model_link(model.model, AutoInvLink()) + @test Effects._model_link(model, AutoInvLink()) == + Effects._model_link(model.model, AutoInvLink()) else invlinks = [iv] end From 5c4c3ac5f80eff4f1e53e35897d579383b9570e8 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 15:05:47 -0500 Subject: [PATCH 13/15] distributions with support outside [0, 1] --- test/delta_method.jl | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/delta_method.jl b/test/delta_method.jl index 3a3e3f7..41ebef4 100644 --- a/test/delta_method.jl +++ b/test/delta_method.jl @@ -40,7 +40,7 @@ using Test @test isapprox(only(eff_emm_trans.upper), 47.5; atol=0.05) @testset "AutoInvLink fails gracefully" begin - # this should work even pre Julia 1.9 because by definition + # this should work even pre Julia 1.9 because by definition # no extension is loaded @test_throws ArgumentError effects(design, model; invlink=AutoInvLink()) end @@ -105,6 +105,20 @@ end end @static if VERSION >= v"1.9" + @testset "Non Link01 GLM link" begin + dat = rdataset("car", "Cowles") + dat[!, :vol] = dat.Volunteer .== "yes" + # this isn't a particularly sensible model, but it's fine for testing + model = glm(@formula(vol ~ Extraversion * Neuroticism), dat, Poisson()) + design = Dict(:Extraversion => [13], + :Neuroticism => [16]) + X = [1.0 13.0 16.0 13 * 16] + eff_manual = effects(design, model; + invlink=Base.Fix1(GLM.linkinv, Link(model.model))) + eff_auto = effects(design, model; invlink=AutoInvLink()) + + @test all(isapprox.(Matrix(eff_manual), Matrix(eff_auto))) + end @testset "link function in a MixedModel" begin model = fit(MixedModel, @formula(use ~ 1 + age + (1 | urban)), From 2dba1a3b7c6eb393365d5259627437d5d922413d Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 16:19:08 -0500 Subject: [PATCH 14/15] phrasing Co-authored-by: Alex Arslan --- ext/EffectsGLMExt.jl | 2 +- src/regressionmodel.jl | 18 ++++++++++++------ 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/ext/EffectsGLMExt.jl b/ext/EffectsGLMExt.jl index 588aa02..24e3844 100644 --- a/ext/EffectsGLMExt.jl +++ b/ext/EffectsGLMExt.jl @@ -3,7 +3,7 @@ module EffectsGLMExt using Effects using GLM: AbstractGLM, Link, Link01, inverselink -using StatsModels: TableRegressionModel +using Effects.StatsModels: TableRegressionModel # TODO: upstream a Link(::TableRegressionModel{<:AbstractGLM}) Effects._model_link(m::TableRegressionModel{<:AbstractGLM}, ::AutoInvLink) = Link(m.model) diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index 0c02298..1f8708f 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -7,9 +7,14 @@ Singleton type indicating that the inverse link should be automatically determined from the model type. -This will only work for model types with an appropriate extension on -Julia 1.9+. If an appropriate extension is not defined, then an error -will occur. +!!! compat "Julia 1.9" + Automatic inverse link determination is implemented using package + extensions, which are available beginning in Julia 1.9. + +An error is thrown if the inverse link cannot be determined. This will +always occur with Julia versions prior to 1.9, and will otherwise occur +when no extension has been loaded that specifies the link function for +the model type. Currently, this is only implemented for GLM.jl and MixedModels.jl """ @@ -71,9 +76,10 @@ Pointwise standard errors are written into the column specified by `err_col`. automatic differentiation. This means that the `invlink` function must be differentiable and should not involve inplace operations. -The special singleton value `AutoInvLink()` can be used to specify that -the appropriate inverse link should be determined and, where possible, a -direct or analytic computation of the derivative is used. +On Julia versions 1.9 or later, the special singleton value `AutoInvLink()` +can be used to specify that the appropriate inverse link should be determined +automatically. In that case, a direct or analytic computation of the derivative +is used when possible. Effects are computed using the model's variance-covariance matrix, which is computed by default using `StatsBas.vcov`. Alternative methods such as the From e8386253cd2001106f75974b46c23a4e5cc5ee10 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 21 Aug 2023 16:19:47 -0500 Subject: [PATCH 15/15] break my tests Co-authored-by: Alex Arslan --- src/regressionmodel.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/regressionmodel.jl b/src/regressionmodel.jl index 1f8708f..91b91bc 100644 --- a/src/regressionmodel.jl +++ b/src/regressionmodel.jl @@ -145,12 +145,13 @@ _invlink_and_deriv(::typeof(identity), η) = (η, 1) # for now, this is private API, but we should see how this goes and whether we can make it public API # so local extensions (instead of Package-Extensions) are better supported _model_link(::RegressionModel, invlink::Function) = invlink -function _model_link(::RegressionModel, ::AutoInvLink) - @static if VERSION < v"1.9" - @error "AutoInvLink requires extensions and is thus not available on Julia < 1.9." +function _model_link(model::RegressionModel, ::AutoInvLink) + msg = string("cannot automatically determine inverse link for models ", + "of type ", typeof(model)) + @static if isdefined(Base, :get_extension) + msg *= "; no appropriate extension has been loaded" end - throw(ArgumentError("No appropriate extension is loaded for automatic " * - "determination of the inverse link for this model type")) + throw(ArgumentError(msg)) end function _difference_method!(eff::Vector{T}, err::Vector{T},