fit_reml

SuperGLM.fit_reml(

X: object,

y: NDArray,

sample_weight: NDArray | None = None,

offset: NDArray | None = None,

*,

max_reml_iter: int = 20,

reml_tol: float | None = None,

pirls_tol: float | None = None,

max_pirls_iter: int | None = None,

lambda2_init: float | None = None,

interaction_mode: str = 'full',

runtime_validation: str | bool = 'auto',

verbose: bool = False,

w_correction_order: int = 1,

) → SuperGLM

Fit with REML estimation of per-term smoothing parameters.

fit_reml() is the smoothness-selection path and does not support a selection penalty: configure selection_penalty=None or 0.0. A Tweedie power within about 0.005 of 1, or further out for repeated responses (whose rows share a lattice phase), makes the dispersion profile multimodal; each of its solves then searches globally (seconds per fit at p = 1.001 on 30,000 distinct rows), and a power too close to 1 for that search raises superglm.NearPoissonDispersionError. It optimizes a Laplace approximate REML objective over per-term smoothing parameters. For sparse/group selection, use fit() or fit_path(). To let REML shrink spline null spaces, use select=True on the spline terms instead.

Parameters:

Xpandas or eager Polars DataFrame

Feature matrix. Lazy frames must be collected before fitting.

yarray-like

Response variable.

sample_weightarray-like, optional

Observation weights, read under the model’s declared weight_semantics. Under "prior" (the default) they state a precision, Var(Y_i | x_i) = phi * V(mu_i) / w_i; under "frequency" they are replication counts and integer weights are likelihood-equivalent to row replication once feature geometry is fixed. The contract sets the likelihood size the REML criterion profiles against, and so reaches the selected smoothing parameters as well as the dispersion. Weights affect fitting but do not enter the linear predictor or automatically scale the conditional mean.

offsetarray-like, optional

Offset term.

max_reml_iterint

Maximum REML outer iterations (default 20).

reml_tolfloat, optional

Stopping tolerance for the smoothing-parameter optimizer. Unset, it resolves per engine: 1e-9 for the Newton engines (exact and discrete), 1e-6 for the step-criterion engines – the SCOP path, and equally the lambda1 selection-penalty path, which optimizes by EFS steps rather than Newton, stops on a per-step lambda-change bound, and does not carry the Newton engines’ determination behavior described below. An explicit value reaches the engine verbatim, except that the discrete engine floors it at 1e-12, below which its cached-W objective differences are numerical noise. Candidate-grade fits (interaction_mode="fast_candidate" and power-search candidates inside estimate_p) pin 1e-6 themselves: they are ranked, never published.

The Newton engines stop when the projected log-lambda gradient and the objective change both fall below reml_tol * (1 + |objective|), judged over the ACTIVE directions: inferentially flat directions are frozen where they stand rather than marched toward the lambda cap. A direction freezes when its gradient falls below max(0.1 * reml_tol, 1e-7) * (1 + |objective|) – a loose tolerance widens that arm; tightening stops at the 1e-7 floor – and its row curvature per penalty dimension (max_j |H_ij| / sqrt(rank_i * rank_j) over the estimated block) is under 1% of the strongest estimated direction’s, anchored at 0.1 per dimension for all-weak models. Rows rather than diagonals so coupled curvature counts; symmetrically per rank so a high-rank random effect and a low-rank spline stay commensurate and a shared cross-term reads the same from both ends; relative to the strongest direction because judging curvature against the objective’s scale would freeze informative directions as the row count grows. The per-direction freeze decision – including the row_curvature, penalty_rank and curvature_bar it judged – is recorded in reml_diagnostics() under profile["reml_freeze_decision"]. Loose values still leave the informative smoothing parameters – and the standard errors computed from them – underdetermined long before predictions are affected: at 1e-6, a converged fit on a synthetic 12k-row stress design publishes a worst-coefficient SE 92% away from the determined answer; 1e-9 pins it to ~0.01%, measured at 3->7 extra Newton iterations on a 97k benchmark fixture and 5->13 on that stress design. A tolerance tighter than the candidate machinery can resolve terminates – on the exact engine, whose line search is where that limit surfaces, and on the discrete engine’s shared-tensor line search – as termination_reason="converged_at_precision" (every active gradient under max(1e-7, reml_tol) * (1 + |objective|), at least one evaluated trial rejected, none left) with converged=True; line_search_failed with converged=False is reserved for genuinely undetermined stalls. The discrete engine’s exits are therefore score_objective_tolerance, active_set_stationary, fixed_lambdas, max_reml_iter and, from a dead line search on a numeric-by-numeric tensor interaction with a known-scale family only (an estimated-scale family takes the generic path even with a tensor interaction), line_search_failed or converged_at_precision, decided by the same predicate the exact engine uses. A dead tensor search whose active gradient is still above that bar is line_search_failed with converged=False; the exit fires only once the candidate’s own working-model step has settled, so the next iteration would repeat this one, and never on the first iteration. On the discrete engine’s other paths a dead search keeps iterating to max_reml_iter: there one working-model update per outer iteration routinely rescues the next search, so a failed search is not evidence of a fixed point. The SCOP engine stops on a per-step lambda-change bound against reml_tol and classifies its endgame as objective_plateau only once steps have stopped contracting – while iterations still buy precision it keeps going toward the strict bound, so an unreachable tolerance ends as an honest plateau classification rather than an early grant.

Versions before 0.20.0 defaulted every engine to 1e-6, so a default-tolerance fit can publish standard errors that differ from 0.19.x by the full determination gap above. The numbers moved once, to the determined values; pass reml_tol=1e-6 to reproduce the old ones.

pirls_tolfloat, optional

Inner PIRLS/IRLS convergence tolerance. Defaults to constructor tol (1e-6). Pass explicitly to override.

max_pirls_iterint, optional

Maximum inner PIRLS iterations per REML step. Defaults to constructor max_iter (100).

lambda2_initfloat, optional

Initial per-group lambda. Defaults to self.lambda2.

interaction_mode{“full”, “fast_candidate”}

"full" runs ordinary REML. "fast_candidate" caps REML outer updates for interaction models and then runs the normal final refit, intended for screening candidate interactions before a full final model fit.

runtime_validation{“auto”, “full”, “skip”} or bool

Controls the post-fit public-runtime parity diagnostic. "auto" validates small fits and skips the full training-row diagnostic for large fits or fast candidate interaction fits. "full" always validates; "skip" skips validation while still canonicalizing the public prediction state.

verbosebool

Print progress.

w_correction_orderint

Order of the W(rho) implicit-differentiation correction. 1 gives the exact objective and gradient with a modified-Newton outer Hessian (default, fast). 2 also includes the available exact d²W/dη² Hessian cross-terms from Wood (2011, Appendix C). Only affects the exact REML path.

Returns:

SuperGLM

The fitted model (self).