-
Notifications
You must be signed in to change notification settings - Fork 4
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Automatic inverse link #56
Changes from all commits
e1eb689
6726f88
ea08c3d
9a447a9
1433b1e
b0355e4
6ac86f8
178f93e
66459e7
b6cf8e2
e79867b
b925065
5c4c3ac
2dba1a3
e838625
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,14 @@ | ||
module EffectsGLMExt | ||
|
||
using Effects | ||
|
||
using GLM: AbstractGLM, Link, Link01, inverselink | ||
using Effects.StatsModels: TableRegressionModel | ||
|
||
# TODO: upstream a Link(::TableRegressionModel{<:AbstractGLM}) | ||
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 |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,9 @@ | ||
module EffectsMixedModelsExt | ||
|
||
using Effects | ||
using MixedModels | ||
using GLM: Link | ||
|
||
Effects._model_link(m::GeneralizedLinearMixedModel, ::AutoInvLink) = Link(m) | ||
|
||
end # module |
Original file line number | Diff line number | Diff line change | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
@@ -2,15 +2,16 @@ using DataFrames | |||||||||||||
using Effects | ||||||||||||||
using GLM | ||||||||||||||
using LinearAlgebra | ||||||||||||||
using MixedModels | ||||||||||||||
using RDatasets: dataset as rdataset | ||||||||||||||
using StableRNGs | ||||||||||||||
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 +38,99 @@ 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 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 | ||||||||||||||
|
||||||||||||||
@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]) | ||||||||||||||
X = [1.0 13.0 16.0 13 * 16] | ||||||||||||||
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 | ||||||||||||||
@testset "invlink = $invlink" for invlink in invlinks | ||||||||||||||
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)) | ||||||||||||||
Comment on lines
+73
to
+75
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
Any reason not to test the approximate equality of the arrays directly? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. because we want element-wise approximate equality, which is a stronger condition (the approximate equality of vectors is based on the norm of the difference vector) |
||||||||||||||
|
||||||||||||||
for level in [0.68, 0.95] | ||||||||||||||
eff = effects(design, model; invlink, level) | ||||||||||||||
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)) | ||||||||||||||
Comment on lines
+83
to
+85
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
Same question as above There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. same answer as above |
||||||||||||||
end | ||||||||||||||
|
||||||||||||||
# 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 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) | ||||||||||||||
|
||||||||||||||
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)) | ||||||||||||||
# 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 | ||||||||||||||
|
||||||||||||||
# 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) | ||||||||||||||
@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()) | ||||||||||||||
|
||||||||||||||
# 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) | ||||||||||||||
@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)), | ||||||||||||||
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 | ||||||||||||||
end | ||||||||||||||
|
||||||||||||||
@testset "identity by another name" begin | ||||||||||||||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
IIUC this will only be loaded if both MixedModels and GLM are loaded. I don't think that restriction is necessary since MixedModels depends on GLM.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
🤷 I wanted to be explicit