Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix lrtest for model families with dispersion #261

Merged
merged 3 commits into from
Aug 5, 2022
Merged

Fix lrtest for model families with dispersion #261

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
094bd2e
Fix `lrtest` for model families with dispersion
nalimilan Aug 2, 2022
ed6f37b
Fixes
nalimilan Aug 5, 2022
655f3ec
Pass `atol`
nalimilan Aug 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 42 additions & 30 deletions src/lrtest.jl
View file
Open in desktop
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
Open in desktop
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