Interaction Screening (PSST)#

SuperGLM.screen_interactions ranks candidate interactions to help choose which models to refit. It probes residual structure using the interaction terms SuperGLM would fit for each pair. Five pair kinds screen today:

ti — spline x spline: a ti()-style tensor deviation surface.
spline_cat — spline x factor: a reference-coded deviation curve for each non-base level.
numeric_cat — numeric x factor: a slope on the numeric for each non-base level.
cat_cat — factor x factor: the cross-level cells of the two-way table.
numeric_numeric — numeric x numeric: a single product tilt.

An OrderedCategorical margin screens as a spline on its mapped level scores — the axis its own refit builds — so an OC x spline pair is a ti row and an OC x factor pair a spline_cat row.

PSST stands for Penalized Smooth Score Test. It uses the fitted mains model to score candidate interactions without fitting a complete new model for each pair. Work includes passes over the observations and matrix calculations in each candidate’s coefficient space; its cost depends on both row count and term size.

model = SuperGLM(
    family="poisson", features=feats, weight_semantics="frequency"
).fit_reml(X, y, sample_weight=exposure)
table = model.screen_interactions(X, y, sample_weight=exposure)
print(table.head())
#  feature_a feature_b kind  statistic      z  edf0  lambda0  n_cells  approx

The test asks, for each pair: after profiling out what the pair’s own main effects already explain, does the working residual carry structure shaped like the block this pair would refit as? For the penalized kinds the probe is evaluated at a ladder of complexity budgets (edf0, default (2, 4, 8, 16) effective degrees of freedom) and each pair is ranked by its best ranking score z, using the candidate’s Gaussian-reference mean and variance. Different budgets probe simple and more flexible shapes. The unpenalized kinds have no penalty to scan: they are evaluated once, at the block’s own dimension.

What the screen does#

Suppose a claims model already has separate smooth effects for driver age and vehicle power. It can express an age effect that applies to every vehicle and a power effect that applies to every driver. An interaction can express something those two effects cannot: for example, high power having a larger effect among young drivers.

PSST asks whether the fitted model’s remaining errors contain that kind of joint structure. For each eligible pair it:

Takes the fitted model’s residual score and working weights.
Builds the interaction shape that SuperGLM would use for that pair, such as a smooth surface or a separate slope for each category.
Removes directions the pair’s own main effects can already represent.
Calculates the candidate’s extra gain in a penalized working quadratic, beyond the pair’s main-effect adjustments. Smooth candidates are tried at several complexity budgets.
Returns a ranking for deciding which complete models to fit and evaluate.

Step 4 solves a small candidate problem while keeping the baseline working model fixed. It does not repeat likelihood fitting and smoothing selection for every pair. The pair projection also does not adjust jointly for every other fitted main effect. This is why the score is a useful probe of a refit, but does not tell us exactly what that refit will achieve.

Thirty features give at most 30 * 29 / 2 = 435 pairs, before unsupported kinds and resource limits are excluded. PSST scores those pairs individually. It does not search all subsets of 435 interactions, assign features to LSS predictors, or find a globally best model. Adding an interaction and refitting can change which pair looks useful next.

How FAST differs#

The original FAST algorithm also ranks interactions against an additive baseline. For a pair of features, it tries one split on each axis, forming four rectangles, and measures how well a constant residual prediction in each rectangle improves the fit. Cumulative tables let it evaluate split locations without scanning all rows again for each split. InterpretML’s current implementation uses binned, objective-specific interaction gains and accepts a baseline through init_score.

PSST instead evaluates the spline, slope or categorical term that the proposed SuperGLM refit would use, with a complexity ladder for penalized terms. That requires more matrix work. The intended benefit is a ranking suited to those refit shapes; superiority for other models or datasets does not follow.

In the recorded comparison, screening ten pairs on 200,000 rows took 2.77–4.89 seconds for PSST and 0.106–0.113 seconds for FAST. These are historical measurements, not timings of every subsequent implementation. Both screens avoid the complete candidate-model refits.

Pair kinds#

kind	probe block	penalty	probe df	refit target
`ti`	centered spline ⊗ centered spline	kron-sum	`edf0` ladder	`TensorInteraction`
`spline_cat`	centered spline menu ⊗ level contrasts	`kron(S_spline, I)`	`edf0` ladder, rank-clamped	`SplineCategorical`
`numeric_cat`	the numeric’s column on each non-base level	none	`L - 1`	`NumericCategorical`
`cat_cat`	level contrasts ⊗ level contrasts	none	`(L1 - 1) * (L2 - 1)`	`CategoricalInteraction`
`numeric_numeric`	the two numerics’ product	none	1	`NumericInteraction`

For the three unpenalized kinds the probe columns are the refit term’s columns: the kron of two contrast menus on a level-pair cell is that cell’s indicator, and the per-level slope columns are exactly what a NumericCategorical builds. For ti and spline_cat the probe spans the refit term’s identifiable deviation space — the part the pair’s own mains cannot absorb.

This holds for cr parents as well as ps, but only because both interaction types now resolve a cubic regression spline the same way. Cubic regression splines carry two bases — a main effect uses the projected-B-spline CubicRegressionSpline, while an interaction marginal uses CardinalCRSpline, which is closer to mgcv’s bs="cr". TensorInteraction had that routing and SplineCategorical did not, so a spline_cat probe could span directions its own refit could not build: on a skewed margin at n_knots=16 the containment of probe span in refit span measured 0.675. Both now go through one shared resolution — including its fallbacks, so a select=True parent or a penalty order the cardinal basis does not implement keeps its own basis on both sides — and the identity holds by construction. See issue #191.

probe df is the unpenalized block’s dimension before the data reduces it: for those kinds the edf0 column reports the rank actually achieved, which is lower when cells are empty or their columns collinear (an 11 x 22 cat_cat table can report 208 rather than 210).

A spline_cat row can also be confirmed as a FactorSmooth when pooling across levels is wanted rather than reference-coded deviations — same parents, penalized level curves; see Interactions.

What gets swept. candidates=None pairs every eligible fitted feature: splines, OrderedCategorical, Categorical and Numeric. Polynomial, RandomEffect and any OrderedCategorical carrying specials= have no screenable margin — a special is a free level with no position on the spline axis, so the pair would need a composite margin. Each of them is reported in table.attrs["deferred_features"], a {feature: reason} mapping, and naming one in candidates raises with that same reason. A Categorical carrying a grouping= is eligible too: screening and the confirmatory interaction refit both validate original labels and apply the same grouping once. Spline x numeric has both margins but no refit target yet, and is deferred until a varying-coefficient term exists.

The deferral has two remedies, and which one applies is a modelling question rather than a screening one. If the Numeric margin is really curved, respec it as a Spline and the pair screens as ti — that is usually the right answer, and the mis-specification was the larger problem anyway. If the margin is genuinely linear, as log1p(Density) is in the example below, respeccing it would be over-fitting to dodge a deferral: leave it, and the pair waits. Reaching for the first remedy by reflex is how a deferral gets “fixed” by introducing a worse mains model. Both cases drop out of the default sweep rather than reporting a null result, and naming one in candidates= raises that deferral specifically, not a generic “unknown feature” error. Pairs the model already fits as an interaction — of any class, FactorSmooth included — are excluded too: the screen profiles only the parent mains, so it cannot re-screen a term already in the model.

A mixed sweep comes back as one sorted table. The example below screens an 80,000-row sample of the freMTPL2 frequency book under weight_semantics="frequency" (y the claim rate ClaimNb / Exposure, sample_weight the exposure, phi estimated at 4.82). Every number here is therefore on the replication contract’s sum(sample_weight) - edf scale, not the default prior one.

The specification is spelled out rather than assumed, because the screen is defined against the fitted mains. A margin mis-specified there does not merely go unscreened — it changes which pairs are screenable at all:

df["LogDensity"] = np.log1p(df["Density"])

features = {
    "DrivAge":    Spline(kind="ps", n_knots=8),
    "VehAge":     Spline(kind="ps", n_knots=12),
    # Bonus-Malus is strongly curved, and its mass piles at the scale
    # minimum; tempered-quantile knots are this library's stated default
    # for exactly that shape (see Feature Types).
    "BonusMalus": Spline(kind="ps", n_knots=12,
                         knot_strategy="quantile_tempered", knot_alpha=0.2),
    # log1p(Density) is the one margin here that genuinely is linear: give
    # it a spline and the smooth collapses to edf 1.4 of rank 11, buying
    # 1.0 deviance. Raw Density is not linear at all -- straightening it
    # is what the log is for.
    "LogDensity": Numeric(),
    "VehBrand":   Categorical(),
    "Region":     Categorical(),
}

table = model.screen_interactions(df, y, sample_weight=exposure)
print(table.to_string(index=False))
#  feature_a  feature_b        kind  statistic         z       edf0      lambda0  n_cells  approx
#     VehAge BonusMalus           ti   5.691019  2.323521   1.999999 1.023168e+02     5244   False
# BonusMalus   VehBrand   spline_cat  11.960557  0.438411   9.999937 7.818059e+09     1012   False
#    DrivAge BonusMalus           ti   1.810771 -0.121587   2.000001 2.611293e+02     7360   False
#     VehAge   VehBrand   spline_cat   7.805321 -0.490671   9.999576 3.779926e+09      627   False
#    DrivAge     VehAge           ti   0.547723 -0.925665   2.000001 5.044218e+01     4560   False
#    DrivAge   VehBrand   spline_cat   1.958228 -1.798192   9.999927 1.301005e+10      880   False
# LogDensity   VehBrand  numeric_cat   1.087366 -1.992925  10.000000 0.000000e+00       11   False
# BonusMalus     Region   spline_cat   6.972190 -2.164533  20.999891 4.198786e+09     2024   False
#    DrivAge     Region   spline_cat   6.753990 -2.198200  20.999860 6.858867e+09     1760   False
# LogDensity     Region  numeric_cat   5.266162 -2.427784  21.000000 0.000000e+00       22   False
#     VehAge     Region   spline_cat   2.008917 -2.930369  20.998740 2.449233e+09     1254   False
#   VehBrand     Region      cat_cat  39.860759 -8.243705 208.000000 0.000000e+00      242   False

Every eligible pair, four kinds, one ranking by z alone. Twelve rows, not fifteen: the three spline x numeric pairs are deferred, so they never enter the sweep.

Read the top two rows together, because they show what z actually ranks. Confirming each by refit — the gate, not the score:

row	kind	`z`	probe df	refit gain
`VehAge x BonusMalus`	`ti`	2.32	2	43.0
`BonusMalus x VehBrand`	`spline_cat`	0.44	10	73.0

The second pair buys more deviance and ranks below the first: 43.0 on 2 df is 21.5 per df against 7.3 for the second pair. The current z adjusts the local score for probe complexity; it does not rank by total refit gain or directly by gain per df. These two gains are measured on the training data. They show that complexity adjustment and total gain can disagree, but do not establish which pair is the better use of a refit budget. That requires held-out gains and refit costs. A holdout study of this ranking, against confirmatory refits on a 200,000-row split, is in Screening Evaluation.

statistic is not comparable down that column: the cat_cat row’s 39.9 is a 208-dimensional block and the numeric_cat row’s 1.1 a 10-dimensional one. The large lambda0 values are bracket edges at clamped rungs, not fitted smoothing parameters. Ten of the twelve rows carry a negative z, which is ordinary: a statistic can land below its rung’s noise floor, and on this book most pairs do. The cat_cat row’s -8.24 is the extreme case — a 208-df block whose statistic, scaled by the whole model’s dispersion, lands at 39.9, so the global phi is conservative for that block. A large negative z only pushes a pair down the queue; it never promotes one. Nothing here was binned or refused (approx is False throughout, no NaN rows).

Read the top row against its own kind’s measured noise maximum below (9.48 for ti): 2.32 does not clear it — and the refit bought 43.0 deviance anyway. That is the screen working as described rather than a contradiction: the floor is the largest value a wide null battery produced, not a threshold a real pair must beat, and the confirmatory refit is what settles the question.

Reading the output#

kind names the block that was probed, and with it the term you would refit to confirm the row. Kinds share one sorted table, but the current z does not give them identical null distributions. Read the reference-law limitation and measured floors below before comparing kinds.
Rank by z, and only z. For the penalized kinds (ti, spline_cat) statistic, edf0 and lambda0 describe the pair’s winning rung, so they are not comparable between rows; at a clamped rung edf0 holds the achieved value and lambda0 is a bracket edge rather than an interpretable smoothing parameter. For the unpenalized kinds (numeric_cat, cat_cat, numeric_numeric) there is a single rung: edf0 reports the block’s achieved rank and lambda0 is 0, and the edf0= argument does not apply to them. statistic is the dispersion-scaled score statistic.
z is a ranking heuristic, not a p-value. By default, the statistic is scaled by the Pearson dispersion of the mains fit; phi= overrides it. Each rung uses the mean and variance of its fixed Gaussian reference law. Estimating the baseline and choosing the best rung affect the sampling distribution, so this standardization does not give a normal score or a p-value. The mathematical distinction and measured floors are below.
n_cells is the grid the probe assembled: the product of the two margins’ grid sizes, where a spline or OC margin contributes its support size, a factor contributes its level count L, and a numeric contributes 1 (so a numeric_cat row reports the factor’s levels and a numeric_numeric row reports 1).
The winning edf0 is a shape diagnostic for penalized rows with material z. A win at rung 2 means tilt-level evidence (a simple in-in surface); wins at 8-16 mean genuinely curved or high-frequency structure. Under the pure null the winning rung is meaningless.
Confirm by refitting. Refit promising pairs as their kind’s refit target and assess held-out predictive gain and cost. Near-tied z values are common. A large score need not imply the greatest predictive gain, and a refit on the same data does not remove selection bias.
NaN rows are skipped or refused pairs, not failures. A gridded pair (ti, spline_cat) is skipped when it exceeds the cell or intermediate budgets even after the quantile-binning fallback, when its tensor curvature block alone is too large for the budget (binning cannot shrink basis dimensions, so those skip with no binning attempted), when its block is too wide to solve inside the budget (see below), or when the statistic degenerates. A numeric_cat pair has no grid to shrink, so a factor too wide for the pair’s blocks is refused rather than approximated: the gate holds the largest of those blocks, the (L + 1)-wide overlap curvature, to the budget — (L + 1)^2 <= max_cells, which admits factors up to 2235 levels — and the block-dimension gate below holds the same width to the solve budget, which at the default is the binding one and admits 1710. Raising max_cells lifts both and computes the pair exactly. numeric_numeric contracts to 3x3 blocks whatever the supports and is never refused. All such rows sort last.
approx=True means the row’s probe basis differs from what a confirmatory refit would build — either a spline margin was quantile-binned for the screen, or that pair’s refit would discretize lossily. Only rows with a spline or OC margin can carry it: cat_cat, numeric_cat and numeric_numeric rows are always approx=False, and refusal is not approximation.

Measured null floors#

These are measured maxima over a null battery, not calibrated quantiles. benchmarks/screening_null_floors.py --seeds 40 fits 160 mains models (four families x 40 seeds, n=8000 rows each) with no interactions in the truth, and screens 3520 pairs. The battery was rerun after correcting the reference variance; all fits completed, with no non-finite scores or warnings. The paired measurement records both normalizations.

kind	rows	mean `z`	p90 `z`	max `z`	probe df
`ti`	480	0.52	1.93	9.48	2-16
`spline_cat`	1440	0.44	1.89	5.66	2-16
`numeric_cat`	960	0.00	1.19	7.54	1-3
`cat_cat`	480	-0.01	1.20	3.98	2-6
`numeric_numeric`	160	-0.07	1.07	4.91	1

Nine of every ten rows fall below 1.1-1.9, depending on kind, yet the maxima reach 3.98-9.48. A z of 5 is therefore not evidence by itself. The correction raises the ti maximum from 7.31 to 9.48 and the spline_cat maximum from 5.53 to 5.66. The three unpenalized kinds retain their normalization. Larger scores after this change do not by themselves demonstrate greater detection power: null scores rise too.

The two largest rows are a rung-2 ti and a 1-df numeric_cat (a slope on a two-level factor). All six largest rows have edf0 <= 4. Low-dimensional quadratics can have substantial right skew even after mean and variance standardization. The unpenalized kinds show this directly:

kind	probe df	rows	max `z`
`numeric_cat`	1	320	7.54
`numeric_cat`	2	320	4.45
`numeric_cat`	3	320	4.78
`cat_cat`	2	160	3.18
`cat_cat`	3	160	3.98
`cat_cat`	6	160	2.90
`numeric_numeric`	1	160	4.91

These maxima are not monotone in df and depend on the number of draws. They do not establish a universal ordering of the kinds’ tails.

The dispersed Gaussian arm has the largest score in four of the five kinds, including both maxima above 6: 7.54 on numeric_cat and 9.48 on ti. The other families reach at most 4.59 (Poisson), 5.34 (gamma) and 5.53 (binomial). A maximum over these four fitted scenarios is not a bound for a new dataset or a different fitted baseline.

Ordered-categorical (OC) and plain spline margins give the following results:

kind	margins	mean `z`	p90 `z`	max `z`
`ti`	plain	0.45	1.90	9.48
`ti`	OC	0.56	1.96	5.57
`spline_cat`	plain	0.48	1.89	5.66
`spline_cat`	OC	0.36	1.89	5.53

The bulk differences are small in this battery, but these measurements do not establish equivalent tails or isolate the effect of a short score grid.

The regression suite uses z < 10 on its specified null fixtures. This is a test bound, not a screening threshold or a probability guarantee. The observed maximum of 9.48 already approaches it, and widening a sweep gives noise more opportunities to produce a large score. Select candidates for refitting using the ranking and the modelling context; assess the refits on held-out data.

What it inherits from the fit#

Screening linearizes at the fitted model: both the offset and the sample_weight used at fit time are applied automatically (weights only when the fit’s were non-unit; pass either only to override), and the mains model’s own smoothness choices define what “leftover” means. Inherited arrays are in training row order, so inheriting requires X/y to be the retained training data — to screen a holdout, subsample, or reordered frame, pass sample_weight (and offset) explicitly. A badly specified mains model screens against the wrong baseline — screening quality is downstream of fit quality. The Pearson dispersion that scales the statistic is attached to the result as table.attrs["phi"] and can be overridden with phi=. By default its residual degrees of freedom follow the fitted model’s declared weight_semantics: sum(sample_weight) - edf under "frequency", the count of positive-weight rows minus edf under "prior" (the default). The screen and fitted model therefore stay on the same dispersion scale. The worked example below declares weight_semantics="frequency", so every number it prints is on the sum(sample_weight) - edf scale.

Factor margins are read through the fitted spec: levels are indexed in the fit’s own order, any LevelGrouping collapse is applied exactly as the fit applies it, and a level the fit never saw raises through the spec’s own validator rather than screening quietly.

Parent smooths must be single-penalty: mains fitted with select=True raise up front, because ti() terms cannot be built on such parents either. For an OrderedCategorical margin the same check applies to its inner spline.

Screening always probes the exact-basis tensor, including for parents that discretize (whose confirmatory refit uses binned marginal supports). That is the same support-discretization gap as the quantile fallback — measured at ~3.5% relative z on signal pairs — and it never affects which basis the confirmatory refit itself uses. To make the gap visible in the output rather than doc-only, any pair whose refit would discretize lossily carries approx=True, applying the gate that refit itself uses: a ti() refit bins its marginal supports only when BOTH parents resolve to fit-time discretization (per-spec discrete overriding the model flag), while a SplineCategorical refit bins whenever its ONE spline parent does. Either way the row is flagged only when some margin that refit would bin has a cardinality exceeding its resolved bin count — lossless binning returns the exact unique support, so low-cardinality rating factors stay approx=False. An OC margin lives on at most one score point per level: it does not reach the binning fallback in practice, and its refits do not discretize at all, so OC pairs stay exact on both sides.

Measured limits#

Corner-localized effects. An interaction confined to a thin corner of the joint support (young driver x high power, with little data there) screens weakly — there is simply little signal-carrying data, and a full refit faces the same limit. Expect low z, not a false positive.
Rare cells and rare levels. A cat_cat cell or a spline_cat level with few rows contributes little to the block, so it screens honestly weak for the same reason — the confirmatory refit is equally starved there. Rare levels cost ranking power, not correctness.
Heavily correlated pairs. At rho ~ 0.85 the joint support is a ridge; the off-ridge tensor directions are thinly identified and a real signal is demoted, not lost. The refit faces the same identifiability.
Continuous x continuous cardinality. Pairs whose unique-value grid (or curvature-intermediate allocation, bounded at a small multiple of the same budget) exceeds max_cells fall back to quantile binning (screen_bins empirical-quantile support points per margin, basis evaluated at within-bin means) and are flagged approx=True in the output. Pairs within budget are always computed exactly — the fallback never touches them. Screening-only: a confirmatory refit of a flagged pair uses the full data.
Block dimension, which costs time rather than memory. max_cells is an allocation ceiling, and every allocation above grows as k^2 in the probe block’s dimension — but the solve grows as k^3, because each rung takes a (k, k) decomposition: a pivoted QR of the pair’s design factor with one SVD beside it. For the gridded kinds a probe column collinear with the overlap span is the routine case rather than the exception (one empty cat_cat cell or one singleton level makes one; the VehBrand x Region row above reports edf0 = 208 against a nominal 210 for exactly this reason), and it is settled by a rank cut on that factor — the Cholesky and its pseudo-inverse fallback this bullet used to describe went with the assembled Grams in issue #257, without moving the cost class. The same knob therefore bounds time: k^3 <= 1000 * max_cells for an unpenalized block, and the same budget against twice the work for a penalized one, whose ladder can bisect rather than clamp. At the default that admits k <= 1709 (cat_cat on two 42-level factors, numeric_cat on a 1710-level factor) and k <= 1357 for ti/spline_cat, measured on the reference box at 0.81 s and 0.67 s per pair — the unpenalized figure a few percent low since that rung’s rank became a count rather than a trace (measured 1.056x at that corner, and it is what makes the reported edf0 the rank rather than k). Those two figures are moment-route figures and the dense path is measured slower since it moved onto design factors (issue #257): 4.09x on the ladder and 4.17x end to end at the widest geometry benchmarked, which carries them to roughly 3.3 s and 2.7 s per pair. The ceilings themselves are dimensional and are unchanged, so no pair is admitted or refused differently; what the move spent is the ~1.5 s per-pair calibration behind them, which is no longer met at the ceiling and is refitted separately. Wider blocks are refused with a NaN row, immediately — binning cannot shrink a basis dimension — and raising max_cells lifts the refusal. For scale, the block the old allocation-only ceiling admitted (k = 4290, two 67-level factors) measured 24 s and 1.3 GB for one pair.
…except spline_cat, which has no block-dimension ceiling. Grouped by level, a spline x categorical pair’s bordered system is an arrow matrix: its curvature is block-diagonal because levels have disjoint row support, its penalty kron(S_spline, I) is block-diagonal on the same grouping, and the only thing coupling the levels is the intercept and the spline main — a border of 1 + k_spline columns, independent of the level count. So a pair the ceiling above refuses is retried through a kernel that factorizes that arrow in time and memory linear in the level count, rather than cubic and quadratic. The dense path still scores every pair it can, unchanged; the structured path only extends where it stops. The block-dimension ceiling is not the only place it stops, so it is not the only door into the kernel: a pair whose support intermediate support * contrasts^2 + L * k_spline^2 blows the same budget is handed over too, whenever the structured path’s transposed intermediate support * k_spline^2 still fits — which happens with the dense block far below the ceiling. Measured: a 22-level factor against a ps(8) margin with 46,119 distinct values is scored by the kernel, exactly, at a block dimension of 231 against the ceiling of 1,357. Exhausting the binning fallback is a third door. So do not read the level count alone to decide which kernel scored a pair — and read the unpinned level count when you do, since a levels= universe with no training rows behind it is pinned to base and widens no block: against a library-default Spline() margin (13 probe columns) over a 400-point support, 140 declared levels of which 105 are populated stays on the dense path at a block of 13 x 104 = 1,352, where 106 populated levels would not. The support matters because the block width is only the first gate — hold that same 140/105 factor and widen the spline and the pair takes the support exit instead at an unchanged block width (measured: 1,800 distinct values dense, 1,848 arrow). Nor is a NaN row a routing signal in reverse: a spline_cat pair whose ladder resolved no direction at any rung gets one on EITHER path. Measured end to end, best of three, one BLAS thread: 200 levels 0.013 s, 1,000 levels 0.048 s, 5,000 levels 0.32 s — against 0.40 s for the 124 levels the dense path tops out at when the block ceiling is what binds. Above roughly five thousand levels the binding cost is no longer screening but fitting the mains model the screen runs against, whose factor contributes one column per level and which no discrete= setting compresses — discretization is support compression for continuous covariates, and a factor is already on its own grid.

Three things bound the structured path, and max_cells scales all three. The level blocks are L x (k_spline + 1)^2, which at the default and a width-11 spline admits 34,722 levels. Before the reference-variance correction, the kernel alone measured there at 1.22 s and 201 MB for a four-rung ladder, and linear below it at 0.16 s for 5,000 and 0.63 s for 20,000. The cell table is support x L, and beside it the spline menu’s outer products are support x k_spline^2; a pair over either is quantile-binned on its spline margin and flagged approx, the same degradation the dense path applies to its own intermediate. And the ladder itself is budgeted in arrow factorizations: two evaluate the bracket, each distinct emitted lambda needs one final variance pass, and a rung inside the bracket needs bisection as well. A ladder clamped to one edge therefore costs three passes. Each final variance pass also adds a QR tree over the levels. A pair that cannot afford the reserved work is refused with a NaN row; the work estimate is not a wall-clock guarantee. Which of those a pair hits depends on its shape — a narrow spline against a huge factor is bounded by the block stacks, a wide spline against a small one by the ladder, since one evaluation is cubic in k_spline where every allocation is quadratic.
A wide factor’s z is an omnibus statistic, and it is diluted. Read the edf0 column before comparing a wide spline_cat row against a narrow one. kron(S_spline, I) leaves the constant and the linear direction unpenalized per level, so its null space is 2(L-1); the mains absorb the constants and roughly L-1 degrees of freedom survive at any penalty. Every rung of the ladder therefore clamps to the same edge, and the pair is only ever tested at edf0 ~ L-1 — the ladder’s whole point, scanning budgets for the one that best matches the signal, is unavailable to it. This is a property of the penalty, not of the level count: it holds for ps, bs and cr margins, whose penalties have a null space, and not for ns, whose penalty is full rank. A spline_cat pair on an ns margin has edf0 = 0 at maximum penalty, so every rung genuinely searches and the ladder costs tens of arrow factorizations. The earlier measurement was 106 against 2 on the same 400-level pair, before final variance passes were added. That difference is why the evaluation budget includes search work. Under a fixed Gaussian reference with d unpenalized identified directions and profiled noncentrality Lambda, the expected score is Lambda / sqrt(2*d). Adding directions while holding that noncentrality fixed reduces the expected score. At fixed dimension and noncentrality, concentrating the interaction in a few levels does not change this omnibus statistic’s distribution. A low score therefore does not establish that every level-specific effect is absent or that a candidate refit would be useless. A pooled or fully penalized refit changes the alternative and its regularization; it does not guarantee stronger screening. Check its benefit against the intended alternative and held-out performance.
Factor and numeric margins have no such cardinality limit. A factor margin never bins — its support is the fitted level set — and a numeric margin never grids at all: it enters its probe linearly, so moments of the numeric accumulated over the other margin’s cells are its exact sufficient statistics, at any cardinality — reduced, level by level, to a factor of the same design rather than to its curvature. Neither margin can raise approx. The one degradation available to a numeric-margin pair is refusal (numeric_cat with a factor wider than its blocks), and a refused row is a NaN row, never an approximated one.

Provenance#

The basic calculation has a direct mathematical interpretation. Let \(U\) be the candidate’s profiled working score, \(V\) its Fisher working curvature and \(S\succeq0\) its penalty. Fix the penalty weight \(\lambda\ge0\). For each candidate coefficient vector \(b\), profile unpenalized adjustments to the intercept and the pair’s main-effect columns. Relative to the nuisance-only profiled optimum, the remaining working gain is

\[ q_\lambda(b)=U^\top b-\tfrac12 b^\top(V+\lambda S)b. \]

When \(V+\lambda S\) is positive definite on the retained candidate space, completing the square gives

\[ \max_b q_\lambda(b)=\tfrac12U^\top(V+\lambda S)^{-1}U=\tfrac12T_\lambda. \]

Thus the raw score is twice the candidate’s extra quadratic gain beyond what those main-effect adjustments can achieve alone. It is not necessarily twice the total gain from the unchanged fitted coefficients. The table reports statistic after division by dispersion, estimated from the mains fit by default or supplied through phi=. The identity is exact for that quadratic; a complete refit changes the working model and need not achieve the same gain.

Classical score testing uses this kind of slope-and-curvature calculation at the null fit, avoiding a fit under every alternative. There is also an established literature on testing smooth components, including Zhang and Lin (2003), and on score-based interaction tests, including GESAT (Lin et al., 2013). These provide related constructions, not a theorem for PSST’s particular penalty choice, pair-only projection and maximum over complexity budgets.

Reference distribution and current limitation#

If the geometry is fixed and \(U\sim N(0,\phi V)\) with known dispersion, the quadratic has a weighted chi-square reference law:

\[ T_\lambda/\phi\overset d=\sum_j a_j Z_j^2,\qquad E(T_\lambda/\phi)=\sum_j a_j=\mathrm{edf}_0,\qquad \operatorname{Var}(T_\lambda/\phi)=2\sum_j a_j^2. \]

Here the \(Z_j\) are independent standard normals and the \(a_j\) are the candidate’s shrinkage eigenvalues. Both execution paths use

\[ z_\lambda=\frac{T_\lambda/\phi-\mathrm{edf}_0}{\sqrt{2\sum_j a_j^2}}. \]

For an identified unpenalized block, every retained \(a_j=1\), so this reduces to the usual \(\sqrt{2\mathrm{edf}_0}\) denominator. With shrinkage, that old denominator was too large. The dense path obtains the sum of squares from its existing decomposition. The structured path computes the same quantity from diagonal and cross-level blocks, without assembling the full smoother.

This standardizes the fixed Gaussian reference at one rung. Fitting the baseline affects the score’s mean and covariance; estimated dispersion and smoothing, non-Gaussian responses, and selecting the best rung introduce further questions. Simulating normal draws from a fixed candidate matrix does not automatically give exact p-values for the fitted-model procedure. PSST currently returns a ranking, with empirical checks of its behavior. Candidate refits and held-out evaluation are the next step.

What PSST combines, and what has been measured#

PSST combines candidate spaces tied to SuperGLM’s refit terms, penalty weights chosen to attain screening-EDF budgets, a maximum over those budgets, and factor-based assembly. Quantized screens are identified by approx. This combination explains the design; establishing its originality requires a more complete literature comparison.

The FAST comparison used the earlier normalization, refits every candidate, and measures held-out gain. In its two specifications, that PSST ranking agrees more closely with held-out refit gain; FAST is faster and agrees more closely with training gain. These are descriptive results from 8 and 10 pairs on one book, not evidence of a universal ranking advantage. Wider candidate sets and replicated data or splits are needed, with the dependence between overlapping pairs accounted for.