"""Common-data model comparisons and explicitly screened likelihood-ratio references.
Nominal parameter counts are not automatically regular model dimensions.
See Self & Liang (1987), JASA 82, 605–610, and Drton (2009), Annals of
Statistics 37, 979–1012, for boundary and singular likelihood-ratio limits.
"""
from __future__ import annotations
import warnings
from collections.abc import Mapping
from dataclasses import dataclass
from typing import Literal
import numpy as np
from scipy.stats import chi2
from knowledgespaces.estimation._data_signature import _data_signature
from knowledgespaces.estimation.blim_em import (
BLIMEstimate,
ConvergenceWarning,
ResponseMatrix,
estimate_blim,
estimate_blim_restarts,
)
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.identifiability import blim_jacobian
from knowledgespaces.estimation.prediction import _positive_integer, predict_blim
from knowledgespaces.estimation.simulate import simulate_blim
from knowledgespaces.estimation.slm import (
SLMEstimate,
estimate_slm,
estimate_slm_restarts,
slm_jacobian,
slm_state_probabilities,
)
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
[docs]
@dataclass(frozen=True)
class ModelScore:
"""Likelihood evaluated on the supplied data; IC penalties are conventional.
``n_parameters`` is the declared free parameter count, not a proven
identifiable dimension. ``BIC`` uses respondent weight; ``BIC_npatterns``
reproduces pks's distinct-observed-pattern convention. Training equality
uses a recorded empirical-measure fingerprint, not only the sample size.
None means an older or manually constructed fit has no fingerprint.
"""
name: str
family: Literal["blim", "slm"]
method: str
n_parameters: int
log_likelihood: float
deviance: float
AIC: float
BIC: float
BIC_npatterns: float
nominal_residual_df: int
pks_residual_df: float
converged: bool
training_data_matches: bool | None
stored_likelihood_matches: bool
[docs]
@dataclass(frozen=True)
class ModelContrast:
"""Adjacent ordered-model contrast; raw signed LR and nominal df difference.
A requested chi-square p-value is NaN when checks fail, with reasons in
``issues``. No reference requested gives None. Numerical ranks at the
null estimate screen local regularity; they do not prove generic/global
identifiability, correct specification, or global likelihood maximization.
"""
null: str
alternative: str
statistic: float
df: int
nested: bool | None
p_value: float | None
issues: tuple[str, ...]
null_rank: int | None = None
alternative_rank_at_null: int | None = None
[docs]
@dataclass(frozen=True)
class ModelComparison:
"""Scores in caller order and contrasts between consecutive models."""
scores: tuple[ModelScore, ...]
contrasts: tuple[ModelContrast, ...]
reference: Literal["none", "chi2"]
n_respondents: float
n_patterns: int
data_signature: str
[docs]
def compare_models(
data: ResponseMatrix,
models: Mapping[str, BLIMEstimate],
*,
reference: Literal["none", "chi2"] = "none",
boundary_tolerance: float = 1e-5,
likelihood_tolerance: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> ModelComparison:
"""Compare BLIM/SLM fits on exactly the supplied complete response data.
Supply at least two distinct names in the desired null-to-alternative
order. Scores are recomputed from parameters and include any held-out or
differently trained model, with training-match flags. MD and MDML fits
may be scored descriptively; their full likelihood is not an ML optimum.
Item columns and state vectors are aligned by labels, not position.
``reference='chi2'`` requests a conditional Wilks reference, screened for
same training measure, converged ML fits, consistent stored likelihoods,
established model nesting, positive nominal dimension difference,
nonnegative LR (up to absolute ``likelihood_tolerance``), interior free
coordinates at the null estimate and full-column-rank response Jacobians.
Free state-support extensions are boundary comparisons and fail this
screen. A failed screen yields NaN plus issues, never a fabricated test.
An accepted screen still needs the scientific assumptions of Wilks's
theorem and adequate optimization/sample size; it is not a proof of them.
Nesting recognizes error equalities/fixed values, fixed/free BLIM state
masses, state-family inclusions and SLM-to-free-prior BLIM restrictions.
General BLIM-to-SLM inclusion is left undetermined. Numerical rank is
evaluated at the null fit in both models, using full response tables;
enumeration/memory guards propagate. MD/MDML and held-out comparisons
remain descriptive even if their numerical scores happen to coincide.
"""
if reference not in ("none", "chi2"):
raise ValueError("reference must be 'none' or 'chi2'.")
if not np.isfinite(boundary_tolerance) or not 0 < boundary_tolerance < 0.5:
raise ValueError("boundary_tolerance must be finite and in (0, .5).")
if not np.isfinite(likelihood_tolerance) or likelihood_tolerance < 0:
raise ValueError("likelihood_tolerance must be finite and nonnegative.")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if len(models) < 2 or any(not isinstance(name, str) or not name for name in models):
raise ValueError("Supply at least two models with nonempty string names.")
data = ResponseMatrix(list(data.items), np.asarray(data.patterns), data.counts)
signature = _data_signature(data)
_, inverse = np.unique(data.patterns, axis=0, return_inverse=True)
counts = np.bincount(inverse, weights=data.effective_counts)
positive = counts[counts > 0]
n, n_patterns = data.n_respondents, len(positive)
saturated = float(np.dot(positive, np.log(positive / n)))
fits = list(models.values())
scores = tuple(
_score(
name,
fit,
data,
signature,
saturated,
n_patterns,
likelihood_tolerance,
max_memory_bytes,
)
for name, fit in models.items()
)
contrasts = []
for index in range(len(fits) - 1):
null, alt = fits[index : index + 2]
left, right = scores[index : index + 2]
nested = _nested(null, alt)
statistic = 2 * (right.log_likelihood - left.log_likelihood)
df = right.n_parameters - left.n_parameters
issues = []
if not np.equal(data.effective_counts, np.floor(data.effective_counts)).all():
issues.append("noninteger_frequencies")
if nested is not True:
issues.append("not_nested" if nested is False else "nesting_undetermined")
if df <= 0:
issues.append("nonpositive_parameter_difference")
if not np.isfinite(statistic):
issues.append("nonfinite_likelihood_ratio")
elif statistic < -likelihood_tolerance:
issues.append("alternative_likelihood_lower")
for score in (left, right):
if score.method != "ML":
issues.append("not_maximum_likelihood")
if not score.converged:
issues.append("unconverged_fit")
if score.training_data_matches is not True:
issues.append(
"training_data_unverified"
if score.training_data_matches is None
else "different_training_data"
)
elif not score.stored_likelihood_matches:
issues.append("stored_likelihood_mismatch")
null_rank = alt_rank = None
p_value = None
if reference == "chi2":
p_value = float("nan")
if not issues:
if not _interior(null, null, boundary_tolerance) or not _interior(
alt, null, boundary_tolerance
):
issues.append("null_on_parameter_boundary")
if not issues:
null_rank = _rank_at(null, null, max_memory_bytes)
alt_rank = _rank_at(alt, null, max_memory_bytes)
if null_rank < left.n_parameters or alt_rank < right.n_parameters:
issues.append("rank_deficient_at_null")
if not issues:
p_value = float(chi2.sf(max(0.0, statistic), df))
contrasts.append(
ModelContrast(
left.name,
right.name,
statistic,
df,
nested,
p_value,
tuple(dict.fromkeys(issues)),
null_rank,
alt_rank,
)
)
return ModelComparison(scores, tuple(contrasts), reference, n, n_patterns, signature)
def _structure(fit: BLIMEstimate) -> KnowledgeStructure:
if type(fit) not in (BLIMEstimate, SLMEstimate):
raise TypeError("Models must be complete-data BLIMEstimate or SLMEstimate objects.")
if len(fit.items) != len(set(fit.items)) or len(fit.states) != len(set(fit.states)):
raise ValueError("Model item and state labels must be unique.")
structure = KnowledgeStructure(fit.items, fit.states)
if structure.states != set(fit.states):
raise ValueError("Model states must already include the empty state and full domain.")
return structure
def _score(
name: str,
fit: BLIMEstimate,
data: ResponseMatrix,
signature: str,
saturated: float,
n_patterns: int,
tolerance: float,
memory: int,
) -> ModelScore:
structure = _structure(fit)
if structure.domain != set(data.items):
raise ValueError("All models must match the response item labels exactly.")
compiled = _compile_constraints(fit.constraints or BLIMConstraints(), fit.items, fit.states)
for values, spec in ((fit.beta, compiled.beta), (fit.eta, compiled.eta)):
if values.shape != (len(fit.items),):
raise ValueError("Model error arrays must match the item labels.")
for group, fixed in zip(spec.groups, spec.fixed, strict=True):
expected = values[group[0]] if fixed is None else fixed
if not np.allclose(values[group], expected, atol=1e-10, rtol=1e-10):
raise ValueError("Fitted error values violate the declared constraints.")
if compiled.pi is not None and not np.allclose(fit.pi, compiled.pi, atol=1e-10, rtol=1e-10):
raise ValueError("Fitted state masses violate the fixed prior.")
if isinstance(fit, SLMEstimate):
if compiled.pi is not None:
raise ValueError("SLM does not support a fixed state prior.")
expected = slm_state_probabilities(structure, fit.g_dict())
if not np.allclose(
expected, [fit.pi_dict()[state] for state in structure], atol=1e-10, rtol=1e-10
):
raise ValueError("SLM state masses do not match its solvability parameters.")
npar = (
compiled.beta.n_free
+ compiled.eta.n_free
+ (
len(fit.items)
if isinstance(fit, SLMEstimate)
else len(fit.states) - 1
if compiled.pi is None
else 0
)
)
ll = 0.0
for start in range(0, data.n_patterns, 1024):
stop = start + 1024
predicted = predict_blim(
structure,
data.patterns[start:stop],
items=data.items,
beta=fit.beta_dict(),
eta=fit.eta_dict(),
pi=fit.pi_dict(),
max_memory_bytes=memory,
)
weights = data.effective_counts[start:stop]
positive = weights > 0
ll += float(np.dot(weights[positive], predicted.log_probabilities[positive]))
training = None if fit.data_signature is None else fit.data_signature == signature
return ModelScore(
name,
"slm" if isinstance(fit, SLMEstimate) else "blim",
fit.method,
npar,
ll,
2 * (saturated - ll),
-2 * ll + 2 * npar,
float(-2 * ll + np.log(data.n_respondents) * npar),
float(-2 * ll + np.log(n_patterns) * npar),
2**data.n_items - 1 - npar,
float(min(2**data.n_items - 1, data.n_respondents) - npar),
fit.converged,
training,
bool(np.isclose(ll, fit.log_likelihood, atol=tolerance, rtol=1e-12)),
)
def _nested(null: BLIMEstimate, alt: BLIMEstimate) -> bool | None:
if set(null.items) != set(alt.items) or not set(null.states) <= set(alt.states):
return False
items = sorted(null.items)
ns = _compile_constraints(null.constraints or BLIMConstraints(), items, null.states)
ats = _compile_constraints(alt.constraints or BLIMConstraints(), items, alt.states)
for n, a in ((ns.beta, ats.beta), (ns.eta, ats.eta)):
group_id = {int(i): k for k, g in enumerate(n.groups) for i in g}
for group, fixed in zip(a.groups, a.fixed, strict=True):
origins = {group_id[int(i)] for i in group}
values = [n.fixed[k] for k in origins]
if fixed is not None:
if any(value != fixed for value in values):
return False
elif len(origins) > 1 and (None in values or len(set(values)) > 1):
return False
if isinstance(alt, SLMEstimate):
return set(null.states) == set(alt.states) if isinstance(null, SLMEstimate) else None
if ats.pi is not None:
if isinstance(null, SLMEstimate) or ns.pi is None:
return False
prior = dict(zip(null.states, ns.pi, strict=True))
return all(
prior.get(state, 0.0) == value for state, value in zip(alt.states, ats.pi, strict=True)
)
return True
def _interior(model: BLIMEstimate, point: BLIMEstimate, tolerance: float) -> bool:
items = sorted(model.items)
spec = _compile_constraints(model.constraints or BLIMConstraints(), items, model.states)
for groups, values in ((spec.beta, point.beta_dict()), (spec.eta, point.eta_dict())):
for group, fixed in zip(groups.groups, groups.fixed, strict=True):
if fixed is None and not tolerance < values[items[int(group[0])]] < 1 - tolerance:
return False
if isinstance(model, SLMEstimate):
assert isinstance(point, SLMEstimate)
return bool(np.all((point.g > tolerance) & (point.g < 1 - tolerance)))
if spec.pi is None:
prior = point.pi_dict()
return all(prior.get(state, 0) > tolerance for state in model.states)
return True
def _rank_at(model: BLIMEstimate, point: BLIMEstimate, memory: int) -> int:
structure = _structure(model)
if isinstance(model, SLMEstimate):
assert isinstance(point, SLMEstimate)
jac = slm_jacobian(
structure,
g=point.g_dict(),
beta=point.beta_dict(),
eta=point.eta_dict(),
constraints=model.constraints,
max_memory_bytes=memory,
)
else:
prior = {state: point.pi_dict().get(state, 0.0) for state in model.states}
jac = blim_jacobian(
structure,
beta=point.beta_dict(),
eta=point.eta_dict(),
pi=prior,
constraints=model.constraints or BLIMConstraints(),
max_memory_bytes=memory,
)
return int(np.linalg.matrix_rank(jac)) if jac.shape[1] else 0
[docs]
@dataclass(frozen=True)
class BootstrapModelComparison:
"""Plug-in null simulation of the specified ML fitting/search procedure.
Observed fits are recomputed with the same settings as every replicate.
Arrays retain all replicates; no unconverged/invalid run is discarded.
``p_value`` is NaN if any selected fit failed convergence or a materially
negative/nonfinite LR occurred. Otherwise it uses (1 + extremes)/(B + 1).
This is a fitted-null Monte Carlo calibration, not an exact finite-sample
test, a proof of bootstrap consistency at singularities, or a guarantee
of global maxima. Both observed and replicate convergence flags are exposed.
"""
null_fit: BLIMEstimate
alternative_fit: BLIMEstimate
statistic: float
statistics: np.ndarray
converged: np.ndarray
p_value: float
n_extreme: int
seed: int | None
n_restarts: int
n_invalid: int
[docs]
def bootstrap_model_comparison(
data: ResponseMatrix,
null_model: BLIMEstimate,
alternative_model: BLIMEstimate,
*,
n_replicates: int = 200,
seed: int | None = None,
n_restarts: int = 1,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> BootstrapModelComparison:
"""Refit two nested model specifications and simulate under the fitted null.
Input fits supply structure, family and constraints; their parameter
values and optimization history are not reused. Both observed models are
fitted anew by ML, then each simulated sample is fitted by the same policy.
One restart uses deterministic standard starts; more use the documented
BLIM/SLM uniform random-start policies. For every alternative fit, an
additional ML run starts at the embedded fitted null, selecting the higher
likelihood (first wins ties). This helps preserve the nesting inequality
without claiming global maximization. Free coordinates can be clipped by
the estimators' documented numerical box at initialization. Starts with
beta + eta >= 1 are allowed, as in the fitted parameter space: an embedded
null is not rejected or reflected to enforce positive discrimination.
Requires established nesting, complete binary data and integer frequency
weights. Every replicate draws N independent respondents from the fitted
null, using its BLIM response probabilities even when it is an SLM. This
can be used for state-support boundary comparisons where ordinary chi2 is
unavailable. It does not settle general nonregular bootstrap validity.
Failed selected fits are retained; an aggregate warning and NaN p-value
require further optimization rather than silently conditioning on success.
"""
n_replicates = _positive_integer(n_replicates, "n_replicates")
n_restarts = _positive_integer(n_restarts, "n_restarts")
max_iter = _positive_integer(max_iter, "max_iter")
memory = _positive_integer(max_memory_bytes, "max_memory_bytes")
if seed is not None and (isinstance(seed, bool) or not isinstance(seed, int) or seed < 0):
raise ValueError("seed must be a nonnegative integer or None.")
if not np.isfinite(tol) or tol <= 0:
raise ValueError("tol must be finite and positive.")
data = ResponseMatrix(list(data.items), np.asarray(data.patterns), data.counts)
for model in (null_model, alternative_model):
if _structure(model).domain != set(data.items):
raise ValueError("Models and data must share item labels.")
if _nested(null_model, alternative_model) is not True:
raise ValueError("Bootstrap comparison requires established null-to-alternative nesting.")
counts = data.effective_counts
if not np.equal(counts, np.floor(counts)).all():
raise ValueError("Bootstrap requires integer frequency counts, not fractional weights.")
n = int(data.n_respondents)
if 8 * n * (6 * data.n_items + 4) + n_replicates * 32 > memory:
raise MemoryError("Bootstrap response generation exceeds max_memory_bytes.")
rng = np.random.default_rng(seed)
def fit_pair(sample: ResponseMatrix) -> tuple[BLIMEstimate, BLIMEstimate]:
null = _refit(null_model, sample, rng, n_restarts, max_iter, tol, memory)
alt = _refit(alternative_model, sample, rng, n_restarts, max_iter, tol, memory)
embedded = _refit(alternative_model, sample, rng, 1, max_iter, tol, memory, start=null)
return null, embedded if embedded.log_likelihood > alt.log_likelihood else alt
with warnings.catch_warnings():
warnings.simplefilter("ignore", ConvergenceWarning)
null_fit, alt_fit = fit_pair(data)
statistic = 2 * (alt_fit.log_likelihood - null_fit.log_likelihood)
statistics = np.empty(n_replicates)
converged = np.empty((n_replicates, 2), dtype=bool)
structure = _structure(null_fit)
for i in range(n_replicates):
sample = simulate_blim(
structure,
n,
beta=null_fit.beta_dict(),
eta=null_fit.eta_dict(),
pi=null_fit.pi_dict(),
seed=rng,
)
assert isinstance(sample, ResponseMatrix)
fitted_null, fitted_alt = fit_pair(sample)
statistics[i] = 2 * (fitted_alt.log_likelihood - fitted_null.log_likelihood)
converged[i] = fitted_null.converged, fitted_alt.converged
valid = np.isfinite(statistics) & (statistics >= -1e-7) & converged.all(axis=1)
observed_valid = (
null_fit.converged and alt_fit.converged and np.isfinite(statistic) and statistic >= -1e-7
)
n_invalid = int((~valid).sum())
n_extreme = int(np.count_nonzero(np.maximum(statistics, 0) >= max(statistic, 0)))
p_value = (
(1 + n_extreme) / (1 + n_replicates) if observed_valid and not n_invalid else float("nan")
)
if not observed_valid or n_invalid:
warnings.warn(
"Model-comparison bootstrap has unconverged or invalid likelihood ratios; "
"all replicates are retained and p_value is NaN. Increase search effort.",
ConvergenceWarning,
stacklevel=2,
)
for array in (statistics, converged):
array.flags.writeable = False
return BootstrapModelComparison(
null_fit,
alt_fit,
statistic,
statistics,
converged,
p_value,
n_extreme,
seed,
n_restarts,
n_invalid,
)
def _refit(
template: BLIMEstimate,
data: ResponseMatrix,
rng: np.random.Generator,
restarts: int,
max_iter: int,
tol: float,
memory: int,
start: BLIMEstimate | None = None,
) -> BLIMEstimate:
structure = _structure(template)
if start is None and restarts > 1:
seed = int(rng.integers(0, 2**32))
if isinstance(template, SLMEstimate):
return estimate_slm_restarts(
structure,
data,
n_restarts=restarts,
seed=seed,
constraints=template.constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=memory,
)
return estimate_blim_restarts(
structure,
data,
n_restarts=restarts,
seed=seed,
constraints=template.constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=memory,
)
beta = 0.1 if start is None else np.array([start.beta_dict()[q] for q in data.items])
eta = 0.1 if start is None else np.array([start.eta_dict()[q] for q in data.items])
if isinstance(template, SLMEstimate):
assert start is None or isinstance(start, SLMEstimate)
return estimate_slm(
structure,
data,
beta_init=beta,
eta_init=eta,
g_init=0.1 if start is None else start.g_dict(),
constraints=template.constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=memory,
)
return estimate_blim(
structure,
data,
beta_init=beta,
eta_init=eta,
pi_init=None
if start is None
else {state: start.pi_dict().get(state, 0.0) for state in structure},
constraints=template.constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=memory,
)