Empirical validation and BLIM diagnostics

Validation against a structure and diagnosis of a fitted probability model answer different questions. Both accept ResponseMatrix, including aggregated counts. Columns align by item labels.

Distances and prerequisite validation

from knowledgespaces.datasets import load_probability
from knowledgespaces.metrics import validate_structure, validate_relation

ds = load_probability(wave="pre", model="K1")
distance = validate_structure(ds.structure, ds.data)
relation = validate_relation(ds.structure, ds.data)
print(distance.di, distance.da)
print(distance.distance_frequencies)  # total weight at distances 0, 1, ...
print(relation.gamma, relation.vc)

For a response pattern \(R\), its distance is \(d(R,\mathcal K)= \min_{K\in\mathcal K}|R\triangle K|\). With weights \(w_R\) and \(N=\sum_R w_R\),

\[ DI=\frac{1}{N}\sum_R w_R d(R,\mathcal K),\qquad DA=\frac{DI}{2^{-|Q|}\sum_{R\subseteq Q} d(R,\mathcal K)}. \]

These implement the definitions documented in kst::kvalidate, citing Schrepp (1999) and Schrepp, Held & Albert (1999). DA uses the uniform response powerset, not a fitted BLIM. Smaller values mean closer data; DA can exceed one and is undefined (nan) for a powerset structure. Neither coefficient is a significance test. The supplied family is used as is, including non-union-closed models such as K1/K2.

pattern_distances follows every input row; zero-weight rows contribute no frequency. uniform_distance_frequencies counts all possible patterns by distance. DA enumerates \(2^{|Q|}\) patterns in chunks. Use compute_da=False for empirical distances and DI without enumeration. Work and memory guards reject oversized requests explicitly.

For every nonreflexive pair (a, b) (b implies a), responses a=1,b=0 contribute to \(N_c\), and a=0,b=1 to \(N_d\). Ties contribute neither.

\[ \gamma=\frac{N_c-N_d}{N_c+N_d},\qquad VC=\frac{N_d}{N\,|S_{\ne}|}. \]

Each respondent can contribute to multiple pairs. Relation generators are transitively closed; equivalent items contribute both directions. Gamma is nan without untied responses; VC is nan without nonreflexive pairs. solution_rates contains fractions in the input item order; multiply by 100 to reproduce kst’s percentages. A general structure’s implied relation describes only pairwise prerequisites, so these indices cannot assess every restriction of that family.

Complete response-table residuals

from knowledgespaces.estimation import BLIMConstraints, estimate_blim

constraints = BLIMConstraints(beta_equal=[ds.data.items], eta_equal=[ds.data.items])
fit = estimate_blim(ds.structure, ds.data, constraints=constraints,
                    max_iter=3000, tol=1e-8)
residuals = fit.residuals(ds.data)
print(residuals.G2, residuals.X2)

from knowledgespaces.viz import plot_blim_residuals  # requires the viz extra
figure = plot_blim_residuals(residuals)
figure.savefig("blim-residuals.png", dpi=150)

For observed count \(O_R\) and expected count \(E_R=N P(R)\), Pearson residuals are \((O_R-E_R)/\sqrt{E_R}\). Signed deviance residuals are

\[ \operatorname{sign}(O_R-E_R) \sqrt{2\{O_R\log(O_R/E_R)-O_R+E_R\}}. \]

The complete table includes unobserved cells and pools duplicate rows. Squared residuals sum to Pearson \(X^2\) and multinomial deviance \(G^2\). This follows the formulas exposed by pks::residuals.blim. A deliberate difference: we retain descriptive residuals when reported degrees of freedom are zero; pks 0.7-0 returns zero deviance residuals.

Cells with \(O=E=0\) contribute zero. An observed, exactly impossible cell gives infinite residuals. Log probabilities distinguish impossibility from underflow. The plot rejects nonfinite residuals explicitly. It displays every cell and a zero reference line; no smoothing or significance threshold is added. It does not reproduce every plotting option in R/MATLAB.

fit.residuals(other_data) also supports held-out responses: their deviance is against the fixed model and can differ from training deviance. Fractional weights support descriptive calculations, without automatically admitting multinomial inference. The complete output table is exponential in domain size even though prediction temporaries are chunked.

Bootstrap parameter replicates

from knowledgespaces.estimation import bootstrap_gof

boot = bootstrap_gof(ds.structure, ds.data, estimate=fit,
                     n_replicates=1000, seed=42, max_iter=3000, tol=1e-8)
summary = boot.parameter_summary()
print(summary.beta_mean, summary.beta_sd)
print(summary.n_used, summary.n_capped)

beta_replicates, eta_replicates and pi_replicates have one row per sample, in the observed fit’s item/state order. g2_replicates, converged_replicates and iterations_replicates retain every refit’s outcome. No unconverged fit is removed by default. parameter_summary(converged_only=True) makes that selection explicit and reports the number used; selection can change the distribution. Standard deviations use ddof=1, giving nan with fewer than two selected replicates. These are empirical dispersions, not automatic confidence intervals. Nonidentifiable models can yield algorithm-dependent parameter dispersion even when predictions agree.

Retaining parameters adds outputs to the existing parametric bootstrap without changing random draws. The authors’ MATLAB toolbox also exposes parameter replicates, but MATLAB was not executed in this comparison and its implementation has not been copied.

Parameter-family trade-offs

from knowledgespaces.estimation import blim_tradeoffs

report = blim_tradeoffs(ds.structure, beta=fit.beta_dict(), eta=fit.eta_dict(),
                       pi=fit.pi_dict(), constraints=constraints)
block = report.block("beta", "eta", "pi")
print(block.rank, block.nullity, block.rank_tolerance)
print(report.parameters)  # named free coordinates, including equality groups

The numerical analysis uses the Jacobian criterion of Stefanutti et al. (2012). Every nonempty beta/eta/pi combination has a JacobianBlock. null_space is an orthonormal basis of directions with numerically zero first derivative of the prediction map. Rows follow report.parameters at block.columns; fixed parameters have no column. The dependent prior mass belongs to the first canonical state.

Each null space includes within-family and cross-family dependencies. Its basis can rotate or change sign between platforms; the span matters. It need not equal the row-reduction basis printed by pks::blimit. The full three-family block detects dependencies absent from pairwise blocks. Rows of the basis show coordinate participation, without identifying a unique set of parameters to fix.

The default SVD cutoff is max(shape) * eps * largest_singular_value for each block; rank_tolerance sets a shared absolute cutoff. This evaluates one point. A deficient derivative at a singular point does not prove generic non-identifiability; at a boundary, null directions need not remain inside the feasible probability region. Use check_identifiability for multiple-point screening and substantive assumptions to justify constraints.

Compare probabilistic models on common data

compare_models evaluates each BLIM/SLM fit on a supplied ResponseMatrix. The default provides descriptive scores and adjacent contrasts in your chosen order; it does not choose a reference distribution automatically.

from knowledgespaces.estimation import compare_models

comparison = compare_models(data, {"restricted": restricted_fit, "full": full_fit})
for score in comparison.scores:
    print(score.name, score.log_likelihood, score.n_parameters, score.AIC, score.BIC)
for contrast in comparison.contrasts:
    print(contrast.statistic, contrast.df, contrast.nested)

Names must be unique nonempty strings and all models must share the response item labels. Columns and state masses align by labels. Scores include full response log likelihood, deviance against the empirical saturated distribution, conventional AIC, BIC using respondent weight, and BIC_npatterns reproducing the pks convention. nominal_residual_df is \(2^{|Q|}-1-p\); pks_residual_df is \(\min(2^{|Q|}-1,N)-p\). Negative values are retained. These are nominal counts, not established identifiable dimensions. A contrast uses \(2(\ell_1-\ell_0)\) and \(p_1-p_0\), never a difference of clipped goodness-of-fit degrees of freedom.

Every new complete BLIM/SLM estimate records data_signature, a fingerprint of its labeled weighted response measure. It ignores row/column order, duplicate-row compression and zero-frequency rows. It detects changed counts that a table-size comparison would miss. It does not identify participants, prove independence, or prove that parameters maximize the likelihood. Differently rounded fractional weights can have different fingerprints. Older/manual estimates may have data_signature=None.

training_data_matches and stored_likelihood_matches distinguish training comparisons from evaluation elsewhere or changed fit objects. Held-out and MD/MDML fits may still be scored descriptively; their scores do not become an ordinary training-data likelihood-ratio test. AIC/BIC penalties for those cases are conventional summaries, not justified automatic selection rules.

Requesting a chi-square reference

Use reference="chi2" to request a Wilks reference. The comparison checks:

  • the same recorded training response measure and integer frequency counts;

  • converged ML estimates, with consistent stored/recomputed likelihoods;

  • established inclusion of the null model in the alternative, and positive nominal parameter difference;

  • a nonnegative likelihood ratio, within the explicit absolute numerical tolerance;

  • interior free coordinates at the null estimate, in both models;

  • full-column-rank response Jacobians at that point in both parameterizations.

If a check fails, p_value is NaN and issues explains why. The default reference="none" has p_value=None and does not compute the boundary/rank screen. The raw signed statistic is retained; negative values within likelihood_tolerance use zero only when evaluating the reference tail. boundary_tolerance is a numerical screen, including the estimators’ clipping limits. Rank checks enumerate the response space and respect memory/domain guards.

Recognized nesting includes fixed/equal error restrictions, fixed/free BLIM priors, inclusion of state families, and SLM restrictions within a BLIM with free state masses. General BLIM-to-SLM inclusion remains undetermined. Adding free-mass states embeds the smaller model at zero probabilities for the added states, so an ordinary chi-square reference is not automatically appropriate. The nominal parameter difference is also insufficient at a rank-deficient point. See Self & Liang (1987) and Drton (2009).

Passing these checks is numerical evidence, not a proof of all Wilks assumptions: adequate optimization, correct specification, independent respondents, prespecified models/constraints and appropriate asymptotics still matter. No post-selection correction or global-identifiability claim is supplied. A constraint estimated from these same data is not a genuinely prespecified fixed null merely because it is entered as beta_fixed or pi_fixed.

Fitted-null bootstrap comparison

from knowledgespaces.estimation import bootstrap_model_comparison

calibration = bootstrap_model_comparison(
    data, restricted_fit, full_fit,
    n_replicates=999, n_restarts=5, max_iter=10000, seed=42,
)
print(calibration.statistic, calibration.p_value, calibration.n_invalid)

The input estimates specify model family, structure and constraints. Both models are fitted anew by ML on the observed data, using the same settings as every simulated sample; template parameter values and their earlier optimization history are not reused. One restart uses standard deterministic starts. More restarts use the documented random-start policies. Every alternative fit also considers a run initialized at the embedded fitted null, then selects the larger likelihood (first on ties). Initializations remain subject to the estimators’ documented numerical box. In particular, a fitted null with beta + eta >= 1 remains an admissible start: its errors are not reflected or altered to impose positive discrimination. Fixed and shared constraints remain in force.

Each replicate draws \(N\) independent complete responses from the fitted null and refits both models. statistics and the two-column converged array retain all runs. If a selected fit is unconverged, or a ratio is nonfinite or materially negative, the p-value is NaN and an aggregate warning requests further search; no failed run is silently discarded. Otherwise the tail is \((1+\#\{T_b\ge T_{obs}\})/(B+1)\), with only numerical negative roundoff set to zero for this tail calculation. n_invalid counts invalid replicates; inspect null_fit and alternative_fit for observed-fit convergence as well.

This is a plug-in calibration under the fitted null and the stated fitting algorithm. It can examine boundary comparisons, but does not prove bootstrap consistency at arbitrary singularities, exact finite-sample validity or global maximization. Models are still required to have established nesting and counts must represent integer frequencies. Incomplete-response model comparisons are not included in this complete-data interface.

cookbook/15_model_comparison.py gives a fully executable regular example: with known zero errors and the two-item powerset, the SLM imposes ordinary 2×2 independence, while the free-prior BLIM is saturated. The LR and its one-degree-of-freedom reference agree with the independent contingency-table calculation; this example does not establish regularity for arbitrary KST models.