LikelihoodRatioTest struct and pretty-printing via coeftable (#265)

* LikelihoodRatioTest struct and pretty-printing via coeftable

* add missing file

* tests

* catch REML in LR test

Author: palday

src/MixedModels.jl

@@ -125,6 +125,7 @@ include("randomeffectsterm.jl")
 include("linearmixedmodel.jl")
 include("gausshermite.jl")
 include("generalizedlinearmixedmodel.jl")
+include("likelihoodratiotest.jl")
 include("linalg/statschol.jl")
 include("linalg/cholUnblocked.jl")
 include("linalg/rankUpdate.jl")

src/likelihoodratiotest.jl

@@ -0,0 +1,106 @@
+"""
+    LikelihoodRatioTest
+
+Results of MixedModels.likelihoodratiotest
+
+## Fields
+* `formulas`: Vector of model formulae
+* `models`: NamedTuple of the `dof` and `deviance` of the models
+* `tests`: NamedTuple of the sequential `dofdiff`, `deviancediff`, and resulting `pvalues`
+
+## Properties
+* `deviance`
+* `pvalues`
+
+"""
+struct LikelihoodRatioTest
+    formulas::AbstractVector{String}
+    models::NamedTuple{(:dof,:deviance)}
+    tests::NamedTuple{(:dofdiff,:deviancediff,:pvalues)}
+end
+
+Base.propertynames(lrt::LikelihoodRatioTest, private = false) = (
+    :deviance,
+    :formulas,
+    :models,
+    :pvalues,
+    :tests,
+)
+
+function Base.getproperty(lrt::LikelihoodRatioTest, s::Symbol)
+    if s == :dof
+        lrt.models.dof
+    elseif s == :deviance
+        lrt.models.deviance
+    elseif s == :pvalues
+        lrt.tests.pvalues
+    elseif s == :formulae
+        lrt.formulas
+    else
+        getfield(lrt, s)
+    end
+end
+
+# backward syntactic but not type compatibility
+Base.getindex(lrt::LikelihoodRatioTest, s::Symbol) = getfield(lrt,s)
+
+"""
+    likelihoodratiotest(m::LinearMixedModel...)
+
+Likeihood ratio test applied to a set of nested models.
+
+Note that nesting of the models is not checked.  It is incumbent on the user to check this.
+"""
+function likelihoodratiotest(m::LinearMixedModel...)
+    if any(getproperty.(getproperty.(m,:optsum),:REML))
+        reduce(==,coefnames.(m))  ||
+                throw(ArgumentError("Likelihood-ratio tests for REML-fitted models are only valid when the fixed-effects specifications are identical"))
+    end
+    m = collect(m)   # change the tuple to an array
+    dofs = dof.(m)
+    formulas = String.(Symbol.(getproperty.(m,:formula)))
+    ord = sortperm(dofs)
+    dofs = dofs[ord]
+    formulas = formulas[ord]
+    devs = objective.(m)[ord]
+    dofdiffs = diff(dofs)
+    devdiffs = .-(diff(devs))
+    pvals = ccdf.(Chisq.(dofdiffs), devdiffs)
+
+    LikelihoodRatioTest(
+        formulas,
+        (dof = dofs, deviance = devs),
+        (dofdiff = dofdiffs, deviancediff = devdiffs, pvalues = pvals)
+    )
+end
+
+function _array_union_nothing(arr::Array{T}) where T
+    Array{Union{T,Nothing}}(arr)
+end
+
+function _prepend_0(arr::Array{T}) where T
+    pushfirst!(copy(arr), -zero(T))
+end
+
+function Base.show(io::IO, lrt::LikelihoodRatioTest; digits=2)
+    println(io, "Model Formulae")
+
+    for (i, f) in enumerate(lrt.formulas)
+        println("$i: $f")
+    end
+    cols = hcat(lrt.models.dof, lrt.models.deviance,
+                _prepend_0(lrt.tests.deviancediff),
+                _prepend_0(lrt.tests.dofdiff),
+                _prepend_0(lrt.tests.pvalues))
+
+    ct = CoefTable(
+        cols, # cols
+        ["model-dof", "deviance", "χ²", "χ²-dof", "P(>χ²)"], # colnms
+        string.(1:length(lrt.formulas)), # rownms
+        5, # pvalcol
+        3 # teststatcol
+    )
+    show(io, ct)
+
+    nothing
+end

src/linearmixedmodel.jl

@@ -487,28 +487,6 @@ end
 
 lowerbd(m::LinearMixedModel) = m.optsum.lowerbd
 
-"""
-    likelihoodratiotest(m::LinearMixedModel...)
-
-Likeihood ratio test applied to a set of nested models.
-
-Note that nesting of the models is not checked.  It is incumbent on the user to check this.
-"""
-function likelihoodratiotest(m::LinearMixedModel...)
-    m = collect(m)   # change the tuple to an array
-    dofs = dof.(m)
-    ord = sortperm(dofs)
-    dofs = dofs[ord]
-    devs = deviance.(m)[ord]
-    dofdiffs = diff(dofs)
-    devdiffs = .-(diff(devs))
-    pvals = ccdf.(Chisq.(dofdiffs), devdiffs)
-    (
-     models = (dof = dofs, deviance = devs),
-     tests = (dofdiff = dofdiffs, deviancediff = devdiffs, p_values = pvals),
-    )
-end
-
 function StatsBase.modelmatrix(m::LinearMixedModel)
     fetrm = first(m.feterms)
     if fetrm.rank == size(fetrm, 2)

test/likelihoodratiotest.jl

@@ -0,0 +1,23 @@
+using MixedModels, Test
+
+@testset "likelihoodratio test" begin
+    slp = MixedModels.dataset(:sleepstudy);
+    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp);
+    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp);
+    lrt = MixedModels.likelihoodratiotest(fm0,fm1);
+
+    @test [deviance(fm0), deviance(fm1)] == lrt.deviance
+    @test deviance(fm0) - deviance(fm1) == first(lrt.tests.deviancediff)
+    @test first(lrt.tests.dofdiff) == 1
+    @test sum(map(length,lrt.tests)) == 3
+    @test sum(map(length,lrt.pvalues)) == 1
+    @test sum(map(length,lrt.models)) == 4
+    @test length(lrt.formulae) == 2
+
+    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp, REML=true);
+    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, REML=true);
+    fm10 = fit(MixedModel,@formula(reaction ~ 1 + days + (1|subj)),slp, REML=true);
+
+    @test_throws ArgumentError MixedModels.likelihoodratiotest(fm0,fm1);
+    lrt =  MixedModels.likelihoodratiotest(fm1,fm10);
+end