Skip to content

Commit 9c68279

Browse files
committed
improve test coverage
1 parent c51c7fc commit 9c68279

File tree

2 files changed

+47
-39
lines changed

2 files changed

+47
-39
lines changed

src/likelihoodratiotest.jl

Lines changed: 40 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -160,26 +160,15 @@ _diff(t::NTuple{N}) where {N} = ntuple(i->t[i+1]-t[i], N-1)
160160
# we only accept one non mixed GLM so that we can require a MixedModel and avoid piracy
161161
function StatsModels.lrtest(m0::GLM_TYPES, m::MixedModel...; atol::Real=0.0)
162162
mods = (m0, m...)
163-
if length(mods) < 2
163+
length(mods) >= 2 ||
164164
throw(ArgumentError("At least two models are needed to perform LR test"))
165-
end
166165
df = dof.(mods)
167-
forward = df[1] <= df[2]
168-
if !all(==(nobs(mods[1])), nobs.(mods))
166+
all(==(nobs(mods[1])), nobs.(mods)) ||
169167
throw(ArgumentError("LR test is only valid for models fitted on the same data, " *
170168
"but number of observations differ"))
171-
end
172-
if forward
173-
for i in 2:length(mods)
174-
if df[i-1] >= df[i] || !isnested(mods[i-1], mods[i]; atol=atol)
175-
throw(ArgumentError("LR test is only valid for nested models"))
176-
end
177-
end
178-
else
179-
for i in 2:length(mods)
180-
if df[i] >= df[i-1] || !isnested(mods[i], mods[i-1], atol=atol)
181-
throw(ArgumentError("LR test is only valid for nested models"))
182-
end
169+
for i in 2:length(mods)
170+
if df[i-1] >= df[i] || !isnested(mods[i-1], mods[i]; atol=atol)
171+
throw(ArgumentError("LR test is only valid for nested models"))
183172
end
184173
end
185174

@@ -238,7 +227,11 @@ function StatsModels.isnested(m1::MixedModel, m2::MixedModel; atol::Real=0.0)
238227
return true
239228
end
240229

241-
function StatsModels.isnested(m1::Union{TableRegressionModel{<:LinearModel},LinearModel}, m2::LinearMixedModel; atol::Real=0.0)
230+
function StatsModels.isnested(m1::TableRegressionModel{Union{<:GeneralizedLinearModel,<:LinearModel}}, m2::MixedModel; atol::Real=0.0)
231+
return _iscomparable(m1.model, m2) && isnested(m1.model, m2; atol)
232+
end
233+
234+
function StatsModels.isnested(m1::LinearModel, m2::LinearMixedModel; atol::Real=0.0)
242235
nobs(m1) == nobs(m2) || return false
243236

244237
size(modelmatrix(m1), 2) <= size(modelmatrix(m2), 2) || return false
@@ -267,15 +260,17 @@ function StatsModels.isnested(m1::GeneralizedLinearModel, m2::GeneralizedLinearM
267260
return true
268261
end
269262

270-
function StatsModels.isnested(m1::TableRegressionModel{<:GeneralizedLinearModel}, m2::GeneralizedLinearMixedModel; atol::Real=0.0)
271-
return isnested(m1.model, m2; atol)
272-
end
273-
274263

275264
#####
276265
##### Helper functions for isnested
277266
#####
278267

268+
"""
269+
_iscomparable(m::LinearMixedModel...)
270+
271+
272+
Check whether LMMs are comparable on the basis of their REML criterion.
273+
"""
279274
function _iscomparable(m::LinearMixedModel...)
280275
isconstant(getproperty.(getproperty.(m, :optsum), :REML)) || throw(
281276
ArgumentError(
@@ -291,47 +286,53 @@ function _iscomparable(m::LinearMixedModel...)
291286
)
292287
end
293288

294-
isconstant(nobs.(m)) ||
295-
throw(ArgumentError("Models must have the same number of observations"))
296-
297289
return true
298290
end
299291

300-
# XXX we need the where clause to distinguish from the general method
301-
# but static analysis complains if we don't use the type parameter
302-
function _samefamily(
303-
::GeneralizedLinearMixedModel{<:AbstractFloat,S}...
304-
) where {S<:Distribution}
305-
return true
306-
end
307-
_samefamily(::GeneralizedLinearMixedModel...) = false
292+
"""
293+
_iscomparable(m::GeneralizedLinearMixedModel...)
308294
295+
296+
Check whethere GLMMs are comparable in terms of their model families and links.
297+
"""
309298
function _iscomparable(m::GeneralizedLinearMixedModel...)
310299
# TODO: test that all models are fit with same fast/nAGQ option?
311300
_samefamily(m...) || throw(ArgumentError("Models must be fit to the same distribution"))
312301

313302
isconstant(string.(Link.(m))) ||
314303
throw(ArgumentError("Models must have the same link function"))
315304

316-
isconstant(nobs.(m)) ||
317-
throw(ArgumentError("Models must have the same number of observations"))
318-
319305
return true
320306
end
321307

308+
"""
309+
_iscomparable(m1::TableRegressionModel, m2::MixedModel)_samefamily(
310+
::GeneralizedLinearMixedModel
311+
312+
Check whethere a TableRegressionModel and a MixedModel have coefficient names indicative of nesting.
313+
"""
322314
function _iscomparable(
323315
m1::TableRegressionModel{<:Union{LinearModel,GeneralizedLinearModel}}, m2::MixedModel
324316
)
325-
_iscomparable(m1.model, m2) || return false
326-
327317
# check that the nested fixef are a subset of the outer
328318
all(in.(coefnames(m1), Ref(coefnames(m2)))) || return false
329319

330320
return true
331321
end
332322

333-
# GLM isn't nested with in LMM and LM isn't nested within GLMM
334-
_iscomparable(m1::Union{LinearModel,GeneralizedLinearModel}, m2::MixedModel) = false
323+
"""
324+
_samefamily(::GeneralizedLinearMixedModel...)
325+
326+
Check whether all GLMMS come from the same model family.
327+
"""
328+
function _samefamily(
329+
::GeneralizedLinearMixedModel{<:AbstractFloat,S}...
330+
) where {S<:Distribution}
331+
# XXX we need the where clause to distinguish from the general method
332+
# but static analysis complains if we don't use the type parameter
333+
return true
334+
end
335+
_samefamily(::GeneralizedLinearMixedModel...) = false
335336

336337

337338
"""

test/likelihoodratiotest.jl

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,7 @@ end
7575
@test !MixedModels.isnested(lm1, fm0)
7676

7777
lrt = likelihoodratiotest(fm0, fm1)
78+
@test lrtest(fm0, fm1) == lrtest(lrt)
7879

7980
@test (deviance(fm0), deviance(fm1)) == lrt.deviance
8081
@test sprint(show, lrt) == "Likelihood-ratio test: 2 models fitted on 180 observations\nModel Formulae\n1: reaction ~ 1 + (1 + days | subj)\n2: reaction ~ 1 + days + (1 + days | subj)\n────────────────────────────────────────────\n DoF -2 logLik χ² χ²-dof P(>χ²)\n────────────────────────────────────────────\n[1] 5 1775.4759 \n[2] 6 1751.9393 23.5365 1 <1e-05\n────────────────────────────────────────────"
@@ -83,6 +84,10 @@ end
8384
lrt = likelihoodratiotest(lm1, fm1)
8485
@test pvalue(lrt) 5.9e-32 atol=1e-16
8586

87+
lrt2 = likelihoodratiotest(lm1.model, fm1)
88+
@test first(lrt2.formulas) == "NA"
89+
@test lrtest(lrt2) == lrtest(lrt)
90+
8691
lrt = likelihoodratiotest(lm0, fm0, fm1)
8792
@suppress @test_throws ArgumentError pvalue(lrt)
8893

@@ -165,6 +170,8 @@ end
165170
)
166171
@suppress @test_throws ArgumentError likelihoodratiotest(gmf, gm_poisson)
167172
@suppress @test_throws ArgumentError likelihoodratiotest(gm0, gm_poisson)
173+
174+
@test !MixedModels._samefamily(gm0, gm_poisson)
168175
# this skips the linear term so that the model matrices
169176
# have the same column rank
170177
@test !MixedModels._isnested(gmf2.mm.m[:, Not(2)], gm0.X)

0 commit comments

Comments
 (0)