LikelihoodRatioTest extends StatsAPI.HypothesisTest (#814)

* LikelihoodRatioTest extends StatsAPI.HypothesisTest

* version bump + news

Author: palday

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v4.33.0 Release Notes
+==============================
+- `LikelihoodRatioTest` now extends `StatsAPI.HypothesisTest` and provides a method for `StatsAPI.pvalue`. [#814]
+
 MixedModels v4.32.0 Release Notes
 ==============================
 - Added `lmm` and `glmm` as convenience wrappers for `fit(LinearMixedModel, ...)` and `fit(GeneralizedLinearMixedModel, ...)` respectively [#810]
@@ -19,7 +23,7 @@ MixedModels v4.29.1 Release Notes
 MixedModels v4.29.0 Release Notes
 ==============================
 - Testbed for experimental support for using PRIMA as an optimization backend introduced via the experimental `prfit!` function. [#799]
-- Julia compat bound raised to current 1.10, i.e. current LTS. [#799] 
+- Julia compat bound raised to current 1.10, i.e. current LTS. [#799]
 
 MixedModels v4.28.0 Release Notes
 ==============================
@@ -612,3 +616,4 @@ Package dependencies
 [#801]: https://github.com/JuliaStats/MixedModels.jl/issues/801
 [#802]: https://github.com/JuliaStats/MixedModels.jl/issues/802
 [#810]: https://github.com/JuliaStats/MixedModels.jl/issues/810
+[#814]: https://github.com/JuliaStats/MixedModels.jl/issues/814

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>"]
-version = "4.32.0"
+version = "4.33.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/MixedModels.jl

@@ -31,7 +31,8 @@ using StaticArrays: StaticArrays, SVector
 using Statistics: Statistics, mean, quantile, std
 using StatsAPI: StatsAPI, aic, aicc, bic, coef, coefnames, coeftable, confint, deviance
 using StatsAPI: dof, dof_residual, fit, fit!, fitted, isfitted, islinear, leverage
-using StatsAPI: loglikelihood, meanresponse, modelmatrix, nobs, predict, r2, residuals
+using StatsAPI:
+    loglikelihood, meanresponse, modelmatrix, nobs, pvalue, predict, r2, residuals
 using StatsAPI: response, responsename, stderror, vcov, weights
 using StatsBase: StatsBase, CoefTable, model_response, summarystats
 using StatsFuns: log2π, normccdf
@@ -129,6 +130,7 @@ export @formula,
     profileσ,
     profilesigma,
     profilevc,
+    pvalue,
     pwrss,
     ranef,
     raneftables,

src/likelihoodratiotest.jl

@@ -15,7 +15,7 @@ Results of MixedModels.likelihoodratiotest
 * `pvalues`
 
 """
-struct LikelihoodRatioTest
+struct LikelihoodRatioTest <: StatsAPI.HypothesisTest
     formulas::AbstractVector{String}
     models::NamedTuple{(:dof, :deviance)}
     tests::NamedTuple{(:dofdiff, :deviancediff, :pvalues)}
@@ -26,6 +26,29 @@ function Base.propertynames(lrt::LikelihoodRatioTest, private::Bool=false)
     return (:deviance, :formulas, :models, :pvalues, :tests)
 end
 
+"""
+    pvalue(lrt::LikelihoodRatioTest)
+
+Extract the p-value associated with a likelihood ratio test.
+
+For `LikelihoodRatioTest`s containing more than one model comparison, i.e. more than two models,
+this throws an error because it is unclear which p-value is desired.
+
+To get p-values for multiple tests, use `lrt.pvalues`.
+"""
+function StatsAPI.pvalue(lrt::LikelihoodRatioTest)
+    pvalues = lrt.pvalues
+    if length(pvalues) > 1
+        throw(
+            ArgumentError(
+                "Cannot extract **only one** p-value from a multiple test result."
+            ),
+        )
+    end
+
+    return only(pvalues)
+end
+
 function Base.getproperty(lrt::LikelihoodRatioTest, s::Symbol)
     if s == :dof
         lrt.models.dof

test/likelihoodratiotest.jl

@@ -58,12 +58,13 @@ end
     @test [deviance(fm0), deviance(fm1)] == lrt.deviance
     @test deviance(fm0) - deviance(fm1) == only(lrt.tests.deviancediff)
     @test only(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 sum(length,lrt.tests) == 3
+    @test sum(length,lrt.pvalues) == 1
+    @test sum(length,lrt.models) == 4
     @test length(lrt.formulae) == 2
     show(IOBuffer(),lrt);
     @test :pvalues in propertynames(lrt)
+    @test only(lrt.pvalues) == pvalue(lrt)
 
     lrt = likelihoodratiotest(lm1,fm1)
     @test lrt.deviance ≈ likelihoodratiotest(lm1.model,fm1).deviance
@@ -73,6 +74,9 @@ end
     @test occursin("-2 logLik", shown)
     @test !occursin("deviance", shown)
 
+    lrt = likelihoodratiotest(lm0, fm0, fm1)
+    @test_throws ArgumentError pvalue(lrt)
+
     # non nested FE between non-mixed and mixed
     @test_throws ArgumentError likelihoodratiotest(lm1, fm0)