diff --git a/src/lrtest.jl b/src/lrtest.jl index 4c616a2d..3e89f6a6 100644 --- a/src/lrtest.jl +++ b/src/lrtest.jl @@ -1,6 +1,7 @@ struct LRTestResult{N} nobs::Int deviance::NTuple{N, Float64} + loglikelihood::NTuple{N, Float64} dof::NTuple{N, Int} pval::NTuple{N, Float64} end @@ -23,9 +24,9 @@ For each sequential pair of statistical models in `mods...`, perform a likelihoo test to determine if the first one fits significantly better than the next. A table is returned containing degrees of freedom (DOF), -difference in DOF from the preceding model, deviance, difference in deviance -from the preceding model, and likelihood ratio and p-value for the comparison -between the two models. +difference in DOF from the preceding model, log-likelihood, deviance, chi-squared +statistic (i.e. absolute value of twice the difference in log-likelihood) +and p-value for the comparison between the two models. Optional keyword argument `atol` controls the numerical tolerance when testing whether the models are nested. @@ -51,23 +52,23 @@ julia> bigmodel = glm(@formula(Result ~ 1 + Treatment + Other), dat, Binomial(), julia> lrtest(nullmodel, model, bigmodel) Likelihood-ratio test: 3 models fitted on 12 observations -────────────────────────────────────────────── - DOF ΔDOF Deviance ΔDeviance p(>Chisq) -────────────────────────────────────────────── -[1] 1 16.3006 -[2] 2 1 15.9559 -0.3447 0.5571 -[3] 4 2 14.0571 -1.8988 0.3870 -────────────────────────────────────────────── +──────────────────────────────────────────────────── + DOF ΔDOF LogLik Deviance Chisq p(>Chisq) +──────────────────────────────────────────────────── +[1] 1 -8.1503 16.3006 +[2] 2 1 -7.9780 15.9559 0.3447 0.5571 +[3] 4 2 -7.0286 14.0571 1.8988 0.3870 +──────────────────────────────────────────────────── julia> lrtest(bigmodel, model, nullmodel) Likelihood-ratio test: 3 models fitted on 12 observations -────────────────────────────────────────────── - DOF ΔDOF Deviance ΔDeviance p(>Chisq) -────────────────────────────────────────────── -[1] 4 14.0571 -[2] 2 -2 15.9559 1.8988 0.3870 -[3] 1 -1 16.3006 0.3447 0.5571 -────────────────────────────────────────────── +──────────────────────────────────────────────────── + DOF ΔDOF LogLik Deviance Chisq p(>Chisq) +──────────────────────────────────────────────────── +[1] 4 -7.0286 14.0571 +[2] 2 -2 -7.9780 15.9559 1.8988 0.3870 +[3] 1 -1 -8.1503 16.3006 0.3447 0.5571 +──────────────────────────────────────────────────── ``` """ function lrtest(mods::StatisticalModel...; atol::Real=0.0) @@ -106,42 +107,53 @@ function lrtest(mods::StatisticalModel...; atol::Real=0.0) end dev = deviance.(mods) - Δdev = _diff(dev) - Δdf = _diff(df) + Δdf = (NaN, _diff(df)...) dfr = Int.(dof_residual.(mods)) - if (forward && any(x -> x > 0, Δdev)) || (!forward && any(x -> x < 0, Δdev)) - throw(ArgumentError("Residual deviance must be strictly lower " * - "in models with more degrees of freedom")) + ll = loglikelihood.(mods) + chisq = (NaN, 2 .* abs.(_diff(ll))...) + + for i in 2:length(ll) + if ((forward && ll[i-1] > ll[i]) || + (!forward && ll[i-1] < ll[i])) && + !isapprox(ll[i-1], ll[i], atol=atol) + throw(ArgumentError("Log-likelihood must not be lower " * + "in models with more degrees of freedom")) + end end - pval = (NaN, chisqccdf.(abs.(Δdf), abs.(Δdev))...) - return LRTestResult(Int(nobs(mods[1])), dev, df, pval) + pval = chisqccdf.(abs.(Δdf), chisq) + return LRTestResult(Int(nobs(mods[1])), dev, ll, df, pval) end function Base.show(io::IO, lrr::LRTestResult{N}) where N Δdf = _diff(lrr.dof) Δdev = _diff(lrr.deviance) + chisq = abs.(2 .* _diff(lrr.loglikelihood)) - nc = 6 + nc = 7 nr = N outrows = Matrix{String}(undef, nr+1, nc) - outrows[1, :] = ["", "DOF", "ΔDOF", "Deviance", "ΔDeviance", "p(>Chisq)"] + outrows[1, :] = ["", "DOF", "ΔDOF", "LogLik", "Deviance", "Chisq", "p(>Chisq)"] outrows[2, :] = ["[1]", @sprintf("%.0d", lrr.dof[1]), " ", - @sprintf("%.4f", lrr.deviance[1]), " ", " "] + @sprintf("%.4f", lrr.loglikelihood[1]), + @sprintf("%.4f", lrr.deviance[1]), + " ", " "] for i in 2:nr outrows[i+1, :] = ["[$i]", @sprintf("%.0d", lrr.dof[i]), @sprintf("%.0d", Δdf[i-1]), - @sprintf("%.4f", lrr.deviance[i]), @sprintf("%.4f", Δdev[i-1]), - string(StatsBase.PValue(lrr.pval[i])) ] + @sprintf("%.4f", lrr.loglikelihood[i]), + @sprintf("%.4f", lrr.deviance[i]), + @sprintf("%.4f", chisq[i-1]), + string(StatsBase.PValue(lrr.pval[i]))] end colwidths = length.(outrows) max_colwidths = [maximum(view(colwidths, :, i)) for i in 1:nc] - totwidth = sum(max_colwidths) + 2*5 + totwidth = sum(max_colwidths) + 2*(nc-1) println(io, "Likelihood-ratio test: $N models fitted on $(lrr.nobs) observations") println(io, '─'^totwidth) diff --git a/test/statsmodel.jl b/test/statsmodel.jl index 37339aab..7e81b6f9 100644 --- a/test/statsmodel.jl +++ b/test/statsmodel.jl @@ -97,6 +97,8 @@ StatsBase.dof_residual(mod::DummyModNoIntercept) = length(mod.y) - length(mod.be StatsBase.nobs(mod::DummyModNoIntercept) = length(mod.y) StatsBase.deviance(mod::DummyModNoIntercept) = sum((response(mod) .- predict(mod)).^2) # isnested not implemented to test fallback +StatsBase.loglikelihood(mod::DummyModNoIntercept) = -sum((response(mod) .- predict(mod)).^2) +StatsBase.loglikelihood(mod::DummyModNoIntercept, ::Colon) = -(response(mod) .- predict(mod)).^2 ## Another dummy model type to test fall-through show method struct DummyModTwo <: RegressionModel @@ -254,15 +256,15 @@ end lr1 = lrtest(m0, m1) @test isnan(lr1.pval[1]) - @test lr1.pval[2] ≈ 0.0010484433450981662 + @test lr1.pval[2] ≈ 3.57538284869704e-6 @test sprint(show, lr1) == """ Likelihood-ratio test: 2 models fitted on 4 observations - ────────────────────────────────────────────── - DOF ΔDOF Deviance ΔDeviance p(>Chisq) - ────────────────────────────────────────────── - [1] 1 14.0000 - [2] 2 1 3.2600 -10.7400 0.0010 - ──────────────────────────────────────────────""" + ────────────────────────────────────────────────────── + DOF ΔDOF LogLik Deviance Chisq p(>Chisq) + ────────────────────────────────────────────────────── + [1] 1 -14.0000 14.0000 + [2] 2 1 -3.2600 3.2600 21.4800 <1e-05 + ──────────────────────────────────────────────────────""" @testset "isnested with TableRegressionModel" begin d = DataFrame(y=y, x1=x1, x2=x2) @@ -275,7 +277,7 @@ end @test StatsModels.isnested(m1, m2) @test StatsModels.isnested(m0, m2) end - + m0 = DummyModNoIntercept(Float64[], ones(4, 0), y) m1 = DummyModNoIntercept([0.3], reshape(x1, :, 1), y) @@ -296,26 +298,36 @@ end "results may not be meaningful"), lrtest(m0, m1)) @test isnan(lr2.pval[1]) - @test lr2.pval[2] ≈ 1.2147224767092312e-5 + @test lr2.pval[2] ≈ 6.128757581368316e-10 # in 1.6, p value printing has changed (JuliaStats/StatsBase.jl#606) if VERSION > v"1.6.0-DEV" @test sprint(show, lr2) == """ Likelihood-ratio test: 2 models fitted on 4 observations - ────────────────────────────────────────────── - DOF ΔDOF Deviance ΔDeviance p(>Chisq) - ────────────────────────────────────────────── - [1] 0 30.0000 - [2] 1 1 10.8600 -19.1400 <1e-04 - ──────────────────────────────────────────────""" + ────────────────────────────────────────────────────── + DOF ΔDOF LogLik Deviance Chisq p(>Chisq) + ────────────────────────────────────────────────────── + [1] 0 -30.0000 30.0000 + [2] 1 1 -10.8600 10.8600 38.2800 <1e-09 + ──────────────────────────────────────────────────────""" else @test sprint(show, lr2) == """ Likelihood-ratio test: 2 models fitted on 4 observations - ────────────────────────────────────────────── - DOF ΔDOF Deviance ΔDeviance p(>Chisq) - ────────────────────────────────────────────── - [1] 0 30.0000 - [2] 1 1 10.8600 -19.1400 <1e-4 - ──────────────────────────────────────────────""" + ────────────────────────────────────────────────────── + DOF ΔDOF LogLik Deviance Chisq p(>Chisq) + ────────────────────────────────────────────────────── + [1] 0 -30.0000 30.0000 + [2] 1 1 -10.8600 10.8600 38.2800 <1e-9 + ──────────────────────────────────────────────────────""" end + + # Test that model with more degrees of freedom that does not improve + # fit compared with simpler model is accepted, even if likelihood is + # lower with some tolerance + lrtest(DummyMod([1], ones(4, 1), y), DummyMod([1, 0], ones(4, 2), y)) + lrtest(DummyMod([1], ones(4, 1), y), DummyMod([1, -1e-8], ones(4, 2), y)) + lrtest(DummyMod([1], ones(4, 1), y), DummyMod([1, -1e-2], ones(4, 2), y), atol=1) + @test_throws ArgumentError lrtest(DummyMod([1], ones(4, 1), y), + DummyMod([1, -1e-2], ones(4, 2), y)) + end