added delta method for glm links

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.5"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -13,6 +14,7 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 DataFrames = "0.22, 1.0, 1.1"
+ForwardDiff = "0.10"
 GLM = "1.5.1"
 MixedModels = "3.9, 4"
 StableRNGs = "1.0"
@@ -26,6 +28,7 @@ julia = "1.6"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
+RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StandardizedPredictors = "5064a6a7-f8c2-40e2-8bdc-797ec6f1ae18"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -34,4 +37,4 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DataFrames", "GLM", "MixedModels", "StableRNGs", "StatsBase", "StatsModels", "StandardizedPredictors", "Test"]
+test = ["DataFrames", "GLM", "MixedModels", "RDatasets", "StableRNGs", "StatsBase", "StatsModels", "StandardizedPredictors", "Test"]

src/Effects.jl

@@ -9,6 +9,7 @@ using StatsBase
 using Tables
 
 using StatsModels: AbstractTerm
+using ForwardDiff
 
 export effects
 export effects!

src/regressionmodel.jl

@@ -1,6 +1,6 @@
 """
     effects!(reference_grid::DataFrame, model::RegressionModel;
-             eff_col=nothing, err_col=:err, typical=mean)
+             eff_col=nothing, err_col=:err, typical=mean, invlink=identity)
 
 Compute the `effects` as specified in `formula`.
 
@@ -25,13 +25,18 @@ If this column does not exist, it is created.
 Pointwise standard errors are written into the column specified by `err_col`.
 
 !!! note
-    Effects are computed on the scale of the transformed response.
-    For models with an explicit transformation, that transformation
-    is the scale of the effects. For models with a link function,
-    the scale of the effects is the _link_ scale, i.e. after
-    application of the link function. For example, effects for
-    logitistic regression models are on the logit and not the probability
-    scale.
+    By default (`invlink=identity`), effects are computed on the scale of the
+    transformed response. For models with an explicit transformation, that
+    transformation is the scale of the effects. For models with a link function,
+    the scale of the effects is the _link_ scale, i.e. after application of the
+    link function. For example, effects for logitistic regression models are on
+    the logit and not the probability scale.
+
+!!! warning
+    If the inverse link is specified via `invlink`, then effects and errors are
+    computed on the original, untransformed scale via the delta method using
+    automatic differentiation. This means that the `invlink` function must be
+    differentiable and should not involve inplace operations.
 
 The reference grid must contain columns for all predictors in the formula.
 (Interactions are computed automatically.) Contrasts must match the contrasts
@@ -54,14 +59,17 @@ Fox, John (2003). Effect Displays in R for Generalised Linear Models.
 Journal of Statistical Software. Vol. 8, No. 15
 """
 function effects!(reference_grid::DataFrame, model::RegressionModel;
-                  eff_col=nothing, err_col=:err, typical=mean)
+                  eff_col=nothing, err_col=:err, typical=mean, invlink=identity)
     # right now this is written for a RegressionModel and implicitly assumes
-    # no link function and the existence of an appropriate formula method
     form = formula(model)
     form_typical = typify(reference_grid, form, modelmatrix(model); typical=typical)
     X = modelcols(form_typical, reference_grid)
     eff = X * coef(model)
     err = sqrt.(diag(X * vcov(model) * X'))
+    if invlink !== identity
+        err .*= ForwardDiff.derivative.(invlink, eff)
+        eff .= invlink.(eff)
+    end
     reference_grid[!, something(eff_col, _responsename(model))] = eff
     reference_grid[!, err_col] = err
     return reference_grid
@@ -80,7 +88,7 @@ end
 """
     effects(design::AbstractDict, model::RegressionModel;
             eff_col=nothing, err_col=:err, typical=mean,
-            lower_col=:lower, upper_col=:upper)
+            lower_col=:lower, upper_col=:upper, invlink=identity)
 
 Compute the `effects` as specified by the `design`.
 
@@ -92,10 +100,11 @@ the lower and upper edge of the error band (i.e. [resp-err, resp+err]).
 """
 function effects(design::AbstractDict, model::RegressionModel;
                  eff_col=nothing, err_col=:err, typical=mean,
-                 lower_col=:lower, upper_col=:upper)
+                 lower_col=:lower, upper_col=:upper, invlink=identity)
+    # TODO: add support for confidence intervals instead of std error
     grid = _reference_grid(design)
     dv = something(eff_col, _responsename(model))
-    effects!(grid, model; eff_col=dv, err_col=err_col, typical=typical)
+    effects!(grid, model; eff_col=dv, err_col, typical, invlink)
     # XXX DataFrames dependency
     grid[!, lower_col] = grid[!, dv] - grid[!, err_col]
     grid[!, upper_col] = grid[!, dv] + grid[!, err_col]

test/delta_method.jl

@@ -0,0 +1,107 @@
+using DataFrames
+using Effects
+using GLM
+using LinearAlgebra
+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])
+    eff_original_scale = effects(design, model; invlink=exp)
+    eff_logscale = effects(design, model)
+    @test all(eff_logscale[!, 3] .≈ log.(eff_original_scale[!, 3]))
+    # the derivative of the exponential function is the exponential function....
+    deriv = exp.(eff_logscale[!, 3])
+    err = eff_logscale.err .* deriv
+    @test all(eff_original_scale.err .≈ err)
+
+    # compare with results from emmeans in R
+    # relatively high tolerances for the point estimates b/c that's very susceptible
+    # to variation in the coef estimates, but the vcov estimates are more stable
+    # emmeans(model,  ~ income * education, level=0.68)
+    eff_emm = effects(Dict(:Income => [6798], :Education => [10.7]), model;
+                      eff_col="log(Prestige)")
+    @test isapprox(only(eff_emm[!, "log(Prestige)"]), 3.84; atol=0.01)
+    @test isapprox(only(eff_emm.err), 0.023; atol=0.005)
+    @test isapprox(only(eff_emm.lower), 3.81; atol=0.005)
+    @test isapprox(only(eff_emm.upper), 3.86; atol=0.005)
+
+    # emmeans(model, ~ income * education, level=0.68, transform="response")
+    eff_emm_trans = effects(Dict(:Extraversion => [12.4], :Neuroticism => [11.5]), model;
+                            invlink=exp, eff_col="Prestige")
+    @test isapprox(only(eff_emm_trans[!, "Prestige"]), 46.4; atol=0.05)
+    @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)
+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])
+    eff = effects(design, model; invlink)
+    X = [1.0 13.0 16.0 13 * 16]
+    # compare with results from GLM.predict
+    pred = DataFrame(predict(model.model, X;
+                             interval=:confidence,
+                             interval_method=:delta,
+                             level=0.68)) # 0.68 is 1 normal quantile, which is just the SE
+    @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)
+    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))
+
+    # 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)
+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)
+    dat = (; :x => x, :y => x .* b1 .+ b0 + randn(StableRNG(42), length(x)))
+    model = lm(@formula(y ~ x), dat)
+
+    invlink = x -> x # same as identity but won't trigger that branch
+    @test invlink !== identity
+
+    design = Dict(:x => 1:20)
+    # if our math on the delta method is correct, then it should work for the
+    # identity case, even though we special case identity() to reduce the
+    # computation
+    eff = effects(design, model; invlink=identity)
+    eff_link = effects(design, model; invlink)
+    # these should be exactly equal b/c derivative is just a bunch of ones
+    # however, we may have to loosen this to approximate equality if
+    # the linear algebra gets very optimized and we start seeing the effects
+    # of associativity in SIMD operations
+    @test eff == eff_link
+end

test/runtests.jl

@@ -7,3 +7,7 @@ end
 @testset "linear regression" begin
     include("linear_regression.jl")
 end
+
+@testset "delta method" begin
+    include("delta_method.jl")
+end