allow specifying CI level for effects() (#50)

* allow effects to take a CI level

* remove done TODO

* fix default in docstring

src/regressionmodel.jl

@@ -144,27 +144,33 @@ end
     effects(design::AbstractDict, model::RegressionModel;
             eff_col=nothing, err_col=:err, typical=mean,
             lower_col=:lower, upper_col=:upper, invlink=identity,
-            vcov=StatsBase.vcov)
+            vcov=StatsBase.vcov, level=nothing)
 
 Compute the `effects` as specified by the `design`.
 
 This is a convenience wrapper for [`effects!`](@ref). Instead of specifying a
 reference grid, a dictionary containing the levels/values of each predictor
 is specified. This is then expanded into a reference grid representing a
 fully-crossed design. Additionally, two extra columns are created representing
-the lower and upper edge of the error band (i.e. [resp-err, resp+err]).
+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.
 """
 function effects(design::AbstractDict, model::RegressionModel;
                  eff_col=nothing, err_col=:err, typical=mean,
                  lower_col=:lower, upper_col=:upper, invlink=identity,
-                 vcov=StatsBase.vcov)
-    # TODO: add support for confidence intervals instead of std error
+                 vcov=StatsBase.vcov, level=nothing)
     grid = expand_grid(design)
     dv = something(eff_col, _responsename(model))
     effects!(grid, model; eff_col=dv, err_col, typical, invlink, vcov)
+    # level=0.68 is approximately one standard error, but it's just enough
+    # off to create all sorts of problems with tests and
+    # cause all our tests to fail, which means it would create problems for
+    # users, so we special case it to maintain legacy behavior
+    level_scale = isnothing(level) ? 1 : sqrt(quantile(Chisq(1), level))
     # XXX DataFrames dependency
-    grid[!, lower_col] = grid[!, dv] - grid[!, err_col]
-    grid[!, upper_col] = grid[!, dv] + grid[!, err_col]
+    grid[!, lower_col] = grid[!, dv] - grid[!, err_col] * level_scale
+    grid[!, upper_col] = grid[!, dv] + grid[!, err_col] * level_scale
     return grid
     # up_low = let dv = getproperty(reference_grid, dv), err = getproperty(reference_grid, err_col)
     #     (; lower_col => dv .- err, upper_col => dv .+ err)

test/delta_method.jl

@@ -46,26 +46,30 @@ end
     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))
+    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))
+
+        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)