Implement HypothesisCoding and deprecate ContrastsCoding (#158)

* add HypothesisCoding, tests, don't export ContrastsCoding

* actually commit tests

* update docs

* add seqdiff coding, update docs and tests

* Apply suggestions from code review

Co-Authored-By: Milan Bouchet-Valat <[email protected]>

* outcome → outcomes

Co-Authored-By: Milan Bouchet-Valat <[email protected]>

* deprecate but export ContrastsCoding; clarify docstring

* fix typo

Co-Authored-By: Milan Bouchet-Valat <[email protected]>

* bump patch version

Co-authored-by: Milan Bouchet-Valat <[email protected]>

Author: kleinschmidt

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsModels"
 uuid = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
-version = "0.6.7"
+version = "0.6.8"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"

docs/src/contrasts.md

@@ -1,5 +1,9 @@
 ```@meta
 CurrentModule = StatsModels
+DocTestSetup = quote
+    using StatsModels
+    using LinearAlgebra
+end
 ```
 
 # Modeling categorical data
@@ -14,7 +18,8 @@ The following contrast coding systems are implemented:
 * [`DummyCoding`](@ref)
 * [`EffectsCoding`](@ref)
 * [`HelmertCoding`](@ref)
-* [`ContrastsCoding`](@ref)
+* [`HypothesisCoding`](@ref)
+* [`SeqDiffCoding`](@ref)
 
 ## How to specify contrast coding
 
@@ -44,13 +49,15 @@ ContrastsMatrix
 DummyCoding
 EffectsCoding
 HelmertCoding
-ContrastsCoding
+SeqDiffCoding
+HypothesisCoding
 ```
 
 ### Special internal contrasts
 
 ```@docs
 FullDummyCoding
+ContrastsCoding
 ```
 
 ## Further details

src/StatsModels.jl

@@ -24,8 +24,10 @@ export
     EffectsCoding,
     DummyCoding,
     HelmertCoding,
+    SeqDiffCoding,
+    HypothesisCoding,
     ContrastsCoding,
-
+    
     coefnames,
     dropterm,
     setcontrasts!,

src/contrasts.jl

@@ -213,7 +213,7 @@ Base.getindex(contrasts::ContrastsMatrix{C,T}, rowinds, colinds) where {C,T} =
 # Making a contrast type T only requires that there be a method for
 # contrasts_matrix(T,  baseind, n) and optionally termnames(T, levels, baseind)
 # The rest is boilerplate.
-for contrastType in [:DummyCoding, :EffectsCoding, :HelmertCoding]
+for contrastType in [:DummyCoding, :EffectsCoding, :HelmertCoding, :SeqDiffCoding]
     @eval begin
         mutable struct $contrastType <: AbstractContrasts
             base::Any
@@ -368,25 +368,153 @@ function contrasts_matrix(C::HelmertCoding, baseind, n)
 end
 
 """
-    ContrastsCoding(mat::Matrix[, base[, levels]])
+    SeqDiffCoding([base[, levels]])
+
+Code each level in order to test "sequential difference" hypotheses, which
+compares each level to the level below it (starting with the second level).
+Specifically, the ``n``th predictor tests the hypothesis that the difference
+between levels ``n`` and ``n+1`` is zero.
+
+# Examples
+
+```jldoctest seqdiff
+julia> seqdiff = StatsModels.ContrastsMatrix(SeqDiffCoding(), ["a", "b", "c", "d"]).matrix
+4×3 Array{Float64,2}:
+ -0.75  -0.5  -0.25
+  0.25  -0.5  -0.25
+  0.25   0.5  -0.25
+  0.25   0.5   0.75
+```
+
+The interpretation of sequential difference coding may be hard to see from the
+contrasts matrix itself.  The corresponding hypothesis matrix shows a clearer
+picture.  From the rows of the hypothesis matrix, we can see that these
+contrasts test the difference between the first and second levels, the second
+and third, and the third and fourth, respectively:
+
+```jldoctest seqdiff
+julia> round.(pinv(seqdiff), digits=2)
+3×4 Array{Float64,2}:
+ -1.0   1.0  -0.0   0.0
+ -0.0  -1.0   1.0  -0.0
+  0.0  -0.0  -1.0   1.0
+```
+
+"""
+SeqDiffCoding
+
+function contrasts_matrix(C::SeqDiffCoding, baseind, n)
+    mat = zeros(n, n-1)
+    for col in 1:n-1
+        mat[1:col, col] .= col-n
+        mat[col+1:end, col] .= col
+    end
+    return mat ./ n
+end
+
+
+"""
+    HypothesisCoding(hypotheses::Matrix[, levels])
+
+Specify how to code a categorical variable in terms of a *hypothesis matrix*.
+For a variable with ``k`` levels, this should be a ``k-1 \times k`` matrix.
+Each row of the matrix corresponds to a hypothesis about the mean
+outcomes under each of the ``k`` levels of the predictor.  The entries in the
+row give the weights assigned to each of these ``k`` means, and the
+corresponding predictor in a regression model estimates the weighted sum of
+these cell means.
+
+For instance, if we have a variable which has four levels A, B, C, and D, and we
+want to test the hypothesis that the difference between the average outcomes for
+levels A and B is different from zero, the corresponding row of the hypothesis
+matrix would be `[-1, 1, 0, 0]`.  Likewise, to test whether the difference
+between B and C is different from zero, the hypothesis vector would be `[0, -1,
+1, 0]`.  To test each "successive difference" hypothesis, the full hypothesis
+matrix would be
+
+```jldoctest hyp
+julia> sdiff_hypothesis = [-1  1  0  0
+                            0 -1  1  0
+                            0  0 -1  1];
+```
+
+Contrasts are derived the hypothesis matrix by taking the pseudoinverse:
+
+```jldoctest hyp
+julia> sdiff_contrasts = pinv(sdiff_hypothesis)
+4×3 Array{Float64,2}:
+ -0.75  -0.5  -0.25
+  0.25  -0.5  -0.25
+  0.25   0.5  -0.25
+  0.25   0.5   0.75
+```
+
+The above matrix is what is produced by constructing a [`ContrastsMatrix`](@ref) from a
+`HypothesisCoding` instance:
+
+```jldoctest hyp
+julia> StatsModels.ContrastsMatrix(HypothesisCoding(sdiff_hypothesis), ["a", "b", "c", "d"]).matrix
+4×3 Array{Float64,2}:
+ -0.75  -0.5  -0.25
+  0.25  -0.5  -0.25
+  0.25   0.5  -0.25
+  0.25   0.5   0.75
+```
+
+The interpretation of the such "sequential difference" contrasts are clear when
+expressed as a hypothesis matrix, but it is not obvious just from looking at the
+contrasts matrix.  For this reason `HypothesisCoding` is preferred for
+specifying custom contrast coding schemes over `ContrastsCoding`.
+
+"""
+mutable struct HypothesisCoding <: AbstractContrasts
+    hypotheses::Matrix
+    contrasts::Matrix
+    base::Nothing
+    levels::Union{Vector,Nothing}
+
+    function HypothesisCoding(hypotheses, base, levels)
+        contrasts = pinv(hypotheses)
+        check_contrasts_size(contrasts, levels)
+        new(hypotheses, contrasts, base, levels)
+    end
+end
+
+HypothesisCoding(mat::Matrix; levels=nothing) =
+    HypothesisCoding(mat, nothing, levels)
+
+function contrasts_matrix(C::HypothesisCoding, baseind, n)
+    check_contrasts_size(C.contrasts, n)
+    C.contrasts
+end
+
+
+"""
+    StatsModels.ContrastsCoding(mat::Matrix[, base[, levels]])
 
 Coding by manual specification of contrasts matrix. For k levels, the contrasts
-must be a k by k-1 Matrix.
+must be a k by k-1 Matrix.  The contrasts in this matrix will be copied directly
+into the model matrix; if you want to specify your contrasts as hypotheses (i.e., 
+weights assigned to each group's cell mean), you should use 
+[`HypothesisCoding`](@ref) instead.
 """
 mutable struct ContrastsCoding <: AbstractContrasts
     mat::Matrix
     base::Any
     levels::Union{Vector,Nothing}
 
     function ContrastsCoding(mat, base, levels)
-        if levels !== nothing
-            check_contrasts_size(mat, length(levels))
-        end
+        Base.depwarn("`ContrastsCoding(contrasts)` is deprecated and will not be exported" *
+                     " in the future, use `HypothesisCoding(pinv(contrasts))` instead.",
+                     :ContrastsCoding)
+        check_contrasts_size(mat, levels)
         new(mat, base, levels)
     end
 end
 
-check_contrasts_size(mat::Matrix, n_lev) =
+check_contrasts_size(mat::Matrix, ::Nothing) = check_contrasts_size(mat, size(mat,1))
+check_contrasts_size(mat::Matrix, levels::Vector) = check_contrasts_size(mat, length(levels))
+check_contrasts_size(mat::Matrix, n_lev::Int) =
     size(mat) == (n_lev, n_lev-1) ||
     throw(ArgumentError("contrasts matrix wrong size for $n_lev levels. " *
                         "Expected $((n_lev, n_lev-1)), got $(size(mat))"))

test/contrasts.jl

@@ -113,6 +113,12 @@
                                         1  1]
     @test coefnames(mf_missing) == ["(Intercept)"; "x: b"]
 
+    # Sequential difference coding
+    setcontrasts!(mf, x = SeqDiffCoding())
+    seqdiff_contr = pinv([-1 1 0
+                          0 -1 1]);
+    @test ModelMatrix(mf).m ≈ hcat(ones(6), seqdiff_contr[[2, 1, 3, 1, 1, 2], :])
+
     # Things that are bad to do:
     # Applying contrasts that only have a subset of data levels:
     @test_throws ArgumentError setcontrasts!(mf, x = EffectsCoding(levels = [:a, :b]))
@@ -125,20 +131,88 @@
     contrasts = [0  1
                  -1 -.5
                  1  -.5]
-    setcontrasts!(mf, x = ContrastsCoding(contrasts))
+    setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts))
     @test ModelMatrix(mf).m == [1 -1 -.5
                                 1  0  1
                                 1  1 -.5
                                 1  0  1
                                 1  0  1
                                 1 -1 -.5]
 
+    contrasts2 = [1 0
+                  1 1
+                  0 1]
+    setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts2))
+    @test ModelMatrix(mf).m == [1  1  1
+                                1  1  0
+                                1  0  1
+                                1  1  0
+                                1  1  0
+                                1  1  1]
+
+    hypotheses2 = pinv(contrasts2)
+    setcontrasts!(mf, x = StatsModels.HypothesisCoding(hypotheses2))
+    @test ModelMatrix(mf).m ≈ [1  1  1
+                               1  1  0
+                               1  0  1
+                               1  1  0
+                               1  1  0
+                               1  1  1]
+
+    # different results for non-orthogonal hypotheses/contrasts:
+    hypotheses3 = [1 1 0
+                   0 1 1]
+    setcontrasts!(mf, x = StatsModels.HypothesisCoding(hypotheses3))
+    @test !(ModelMatrix(mf).m ≈ [1  1  1
+                                 1  1  0
+                                 1  0  1
+                                 1  1  0
+                                 1  1  0
+                                 1  1  1])
+
     # throw argument error if number of levels mismatches
-    @test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(contrasts[1:2, :]))
-    @test_throws ArgumentError setcontrasts!(mf, x = ContrastsCoding(hcat(contrasts, contrasts)))
+    @test_throws ArgumentError setcontrasts!(mf, x = StatsModels.ContrastsCoding(contrasts[1:2, :]))
+    @test_throws ArgumentError setcontrasts!(mf, x = StatsModels.ContrastsCoding(hcat(contrasts, contrasts)))
 
     # contrasts types must be instantiated (should throw ArgumentError, currently
     # MethodError on apply_schema)
     @test_broken setcontrasts!(mf, x = DummyCoding)
 
+    @testset "hypothesis coding" begin
+
+        # to get the scaling right, divide by three (because three levels)
+        effects_hyp = [-1 2 -1
+                       -1 -1 2] ./ 3
+
+        @test modelmatrix(setcontrasts!(mf, x = HypothesisCoding(effects_hyp))) ≈
+            modelmatrix(setcontrasts!(mf, x = EffectsCoding()))
+
+        d2 = DataFrame(y = rand(100),
+                       x = repeat([:a, :b, :c, :d], inner=25))
+
+        sdiff_hyp = HypothesisCoding([-1 1 0 0
+                                      0 -1 1 0
+                                      0 0 -1 1])
+
+        effects_hyp = HypothesisCoding([-1 3 -1 -1
+                                        -1 -1 3 -1
+                                        -1 -1 -1 3] ./ 4)
+
+        f = apply_schema(@formula(y ~ 1 + x), schema(d2))
+
+        f_sdiff = apply_schema(f, schema(d2, Dict(:x => sdiff_hyp)))
+        f_effects = apply_schema(f, schema(d2, Dict(:x => effects_hyp)))
+
+        y_means = by(d2, :x, :y => mean).y_mean
+        
+        y, X_sdiff = modelcols(f_sdiff, d2)
+        @test X_sdiff \ y ≈ [mean(y_means); diff(y_means)]
+
+        y, X_effects = modelcols(f_effects, d2)
+        @test X_effects \ y ≈ [mean(y_means); y_means[2:end] .- mean(y_means)]
+
+        @test X_effects ≈ modelcols(apply_schema(f.rhs,
+                                                 schema(d2, Dict(:x=>EffectsCoding()))),
+                                    d2)
+    end
 end