Skip to content

Commit

Permalink
Fix lrtest for model families with dispersion (#261)
Browse files Browse the repository at this point in the history 
`lrtest` relied on the deviance rather than the log-likelihood, which is
not correct for model families where a dispersion parameter needs to be
taken into account. Scaling the deviance would be more efficient than
computing the log-likelihood, but there is currently no generic API for this
and this may not work for non-GLM models, so simply call `loglikelihood`.
We could imagine defining a
`likelihoodratio(m1, m2) = loglikelihood(m1) - loglikelihood(m2)` method
that packages could override for performance, but this may not be worth it.

Also make the check that more complex models have a strictly better fit than
simpler nested ones less strict. The more complex model may have the same
deviance, and due to approximations it may even have a slightly higher deviance.
  • Loading branch information
@nalimilan
nalimilan authored Aug 5, 2022
1 parent 5299399 commit 24a4e47
Show file tree
Hide file tree
Showing 2 changed files with 75 additions and 51 deletions.
72 changes: 42 additions & 30 deletions src/lrtest.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
54 changes: 33 additions & 21 deletions test/statsmodel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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

0 comments on commit 24a4e47

Please sign in to comment.