mixle.inference.glm module

Generalized linear models and penalized / robust / quantile regression on plain arrays.

A array-level regression toolkit (operating on a design matrix X and response y, independent of the PPL DSL in mixle.ppl.regression):

  • glm() – exponential-family GLMs by iteratively reweighted least squares, with explicit family/link objects (Gaussian, Binomial, Poisson, Gamma, inverse-Gaussian, negative-binomial), offsets, prior weights under a stated convention (weight_type: analytic precision weights by default, frequency counts on request), and optional sandwich (robust) standard errors in the HC0-HC3 leverage-corrected variants.

  • ridge_regression(), elastic_net() (and lasso()) – L2 / L1 / mixed penalised linear regression; the elastic net is solved by coordinate descent.

  • robust_regression() – Huber / Tukey M-estimation, down-weighting outliers via IRLS on a robust scale.

  • quantile_regression() – the conditional tau-quantile by IRLS on the check loss.

glm returns a result with coefficients, standard errors, and Wald inference; the penalized, robust, and quantile fits return coefficients and fits WITHOUT standard errors (their correct inferential machinery – debiased lasso, M-estimator sandwiches, quantile sparsity estimates – is deliberately not faked here; audit G-6).

class Link(name, g, inv, mu_eta)[source]

Bases: object

A link function eta = g(mu) with its inverse and derivative dmu/deta.

Parameters:
class Family(name, variance, canonical, unit_deviance, estimate_dispersion, extra=1.0, default_link=None)[source]

Bases: object

An exponential-family error model: variance function, canonical link, deviance, dispersion.

canonical names the family’s mathematical canonical link; default_link (when set) is the link glm() fits with when none is requested. They are separate so families whose canonical link is numerically awkward (gamma’s inverse, inverse-Gaussian’s inverse_squared – neither keeps mu positive) can default to log without mislabeling the canonical link.

Parameters:
class GLMResult(coef, se, fitted, deviance, dispersion, log_likelihood, n_iter, family, link, cov, converged=True, rank=None, residual_df=None, dispersion_estimated=False, weight_type='analytic', ll_nobs=None, separated=False, cov_type='model', _link=None)[source]

Bases: object

Fitted GLM.

Parameters:
coef

(p,) coefficient estimates.

Type:

numpy.ndarray

se

(p,) standard errors (model-based, or the robust variant named by cov_type). All NaN on a saturated dispersion-estimating fit – see dispersion.

Type:

numpy.ndarray

cov_type

which estimator se / cov came from: "model" for phi (X'WX)^-1, or "HC0" | "HC1" | "HC2" | "HC3" for the sandwich variants (see glm()’s robust argument). Recorded on the result for the same reason weight_type is: the meaning of a headline output must be readable off the result.

Type:

str

fitted

(n,) fitted means mu.

Type:

numpy.ndarray

deviance

residual deviance.

Type:

float

dispersion

estimated/assumed dispersion phi – the residual_df-corrected (Pearson) estimate the standard errors use, NOT the value the log-likelihood is evaluated at. NaN when a dispersion-estimating family is fit with residual_df <= 0 (a SATURATED fit – as many estimated parameters as observations): the coefficients, fitted values and deviance are still returned, but the dispersion, every standard error and the likelihood are undefined there, announced by a UserWarning at fit time, and z_values / p_values / aic / bic refuse with that reason.

Type:

float

residual_df

residual degrees of freedom, n - rank, where n counts rows under weight_type="analytic" and sum(weights) under weight_type="frequency". This is the divisor of dispersion and the t reference p_values uses, so it is where the weight convention shows up in the standard errors.

Type:

int | None

log_likelihood

maximised log-likelihood, or None when the supplied family/response does not define one – or when the fit is saturated (see dispersion), where the likelihood is unbounded in the dispersion and so has no maximised value. For dispersion-estimating families this is evaluated at the family’s own dispersion MLE (see _dispersion_mle), which differs from dispersion by design: the covariance wants the df-corrected estimate, the likelihood wants its maximiser. It is also computed under weight_type: analytic weights enter as per-observation precisions (Var = phi / w_i), frequency weights as replicate counts.

Type:

float | None

ll_nobs

how many observations log_likelihood is a sum over – the sample size bic penalises with. Row count for an analytic-weighted dispersion-estimating family; sum(weights) whenever the likelihood is a weighted sum of per-unit densities (every frequency-weighted fit, and every fixed-dispersion family, whose discrete densities have no analytic-weight form). Pairing a likelihood over sum(w) units with a log(rows) penalty is what made BIC agree with neither convention.

Type:

float | None

weight_type

the convention weights were read under; see glm().

Type:

str

separated

True when the binomial design separates the classes, so at least one coefficient has no finite maximum-likelihood estimate. deviance, fitted and log_likelihood are still correct (compare nested models by deviance), but coef and se along the separating direction record only where IRLS stopped, so z_values / p_values refuse. Announced by a UserWarning at fit time.

Type:

bool

n_iter

IRLS iterations to convergence.

Type:

int

converged

whether the IRLS convergence criterion was met. Public fits currently raise instead of returning this as false.

Type:

bool

rank

effective NUMERICAL rank of the weighted design sqrt(W) X at the solution, from the same SVD the covariance is computed from (cutoff at cond(X), not cond(X)^2; audit B4). rank < p – exact or near collinearity – is announced by a UserWarning at fit time and makes z_values / p_values refuse.

Type:

int | None

family / link

names.

cov

(p, p) coefficient covariance.

Type:

numpy.ndarray

Call forms: aic / bic are PROPERTIES (fit.aic, no parentheses) while z_values / p_values are METHODS (fit.z_values()). Both call forms are frozen for compatibility – changing either direction breaks every existing caller – so the split is stated here instead: mixing them up fails with TypeError: 'float' object is not callable (calling the property) or a bound-method object where an array was expected (reading the method), and neither traceback names the actual mistake. All four refuse with a stated reason (no likelihood, rank deficiency, separation, saturation) rather than return NaN.

predict(x, *, offset=None)[source]

Predict the mean response mu at new design rows x.

Parameters:
Return type:

ndarray

property aic: float

-2 ll + 2 k with k = rank + estimated dispersion.

Type:

Akaike information criterion

property bic: float

Bayesian information criterion over the observations the log-likelihood spans.

The penalty has to count the SAME observations the likelihood does (ll_nobs). Deriving it from residual_df instead pairs a likelihood summed over sum(weights) replicates with a log(rows) penalty on a frequency table, which is correct under neither weight convention – 107.7 where the two defensible answers were 116.5 and a likelihood of -3.5.

z_values()[source]

Return Wald statistics for fitted coefficients (full-rank designs only).

Return type:

ndarray

p_values()[source]

Two-sided Wald p-values: Student-t on the residual df when the dispersion was ESTIMATED (Gaussian/Gamma/inverse-Gaussian – the plug-in-dispersion normal reference rejected 9% at nominal 5% at n=8; audit G-1), normal otherwise (fixed-dispersion families), both asymptotic in the non-Gaussian mean model.

Return type:

ndarray

exception PerfectSeparationError[source]

Bases: RuntimeError

The binomial classes are perfectly separated (or perfectly predicted): no finite MLE exists.

Raised by glm() when IRLS runs the fitted binomial probabilities all the way onto 0 or 1 and can no longer take a step – the terminal form of (quasi-)complete separation, where |coef| diverges and the working weights vanish. A subclass of RuntimeError, so pre-existing handlers keep working.

Separation that stops the iteration just SHORT of that point is the same statistical condition and is reported too, but as a UserWarning plus GLMResult.separated, because there the fit still carries a correct deviance / fitted / likelihood; only its Wald branch is undefined and refuses. Which of the two a given data set produces is a numerical accident, so never read the absence of this exception as evidence of no separation – read separated.

glm(x, y, *, family='gaussian', link=None, offset=None, weights=None, weight_type='analytic', max_iter=100, tol=1e-8, robust=False, theta=1.0)[source]

Fit a generalized linear model by iteratively reweighted least squares.

Parameters:

  • x (ndarray) – (n, p) design matrix (include an intercept column explicitly if wanted).

  • y (ndarray) – (n,) response (counts, 0/1 or proportions, positive reals, … per the family).

  • family (str | Family) – "gaussian", "binomial", "poisson", "gamma", "inverse_gaussian", "negativebinomial", or a Family.

  • link (str | Link | None) – a link name, or a Link instance to use directly (e.g. a custom link); defaults to the family’s default link (the canonical link, except gamma / inverse-Gaussian which default to log – their canonical inverse / inverse_squared links do not keep mu positive).

  • offset (ndarray | None) – (n,) known additive term on the linear-predictor scale (e.g. log exposure).

  • weights (ndarray | None) – (n,) non-negative prior weights, at least one positive. Zero-weight rows drop out of the fit entirely. What a weight MEANS is set by weight_type – the two readings agree on the coefficients and on nothing else.

  • weight_type (str) –

    the convention weights are read under.

    • "analytic" (default; also called prior, precision or variance weights – the McCullagh–Nelder reading, and R’s for glm): row i is ONE observation measured with precision w_i, so Var(y_i) = phi V(mu_i) / w_i. The sample size is the number of rows: residual_df = rows - rank, dispersion divides by that, and the log-likelihood is that of n independent observations with dispersions phi/w_i. Where the dispersion is ESTIMATED, phi then absorbs the weights’ scale, so multiplying every weight by a constant leaves the whole fit unchanged – coefficients, standard errors, log-likelihood, AIC and BIC alike – and only the relative precisions matter. (Not so where phi is fixed at 1: for binomial / Poisson / negative-binomial the weight scale IS the amount of information, and halving every weight widens the standard errors by sqrt(2).) Use this for inverse-variance weighting, measurement precisions or sampling weights; fractional weights are accepted.

    • "frequency": w_i is a COUNT – row i stands for w_i identical observations – so it must be a whole number. Every reported quantity then equals what fitting the expanded sum(w)-row data set gives: residual_df = sum(w) - rank, dispersion divides by that, robust standard errors replicate the score, and AIC/BIC span sum(w) observations. Use it for aggregated frequency tables – survey cells, contingency tables, “distinct rows plus a count” extracts.

    The gap is not cosmetic. On a 6-row table whose counts sum to 114, the same call gives se = [0.369, 0.135] analytically and [0.070, 0.025] as frequencies – a factor sqrt((114 - 2) / (6 - 2)) = 5.29 – with p-values of 5e-04 against 6e-82. Fixed dispersion families (binomial, Poisson, negative-binomial) return identical coefficients AND standard errors either way, because phi is never estimated there; only residual_df and the BIC sample size move. Their likelihood is sum w log f under both conventions (a discrete density has no analytic-weight form), so ll_nobs is sum(w) for them regardless – see GLMResult.ll_nobs.

  • max_iter (int) – IRLS controls (convergence on the relative deviance change).

  • tol (float) – IRLS controls (convergence on the relative deviance change).

  • robust (bool | str) –

    which coefficient-covariance estimator to report (recorded on the result as cov_type). False – the default – is the model-based phi (X'WX)^-1; True or "HC0" is the Huber–White sandwich; "HC1" | "HC2" | "HC3" apply the standard finite-sample leverage corrections to the same PSD-by-construction sandwich (MacKinnon–White 1985):

    • HC0: the plain sandwich. No finite-sample correction, so on small or leveraged designs it is biased LOW – measured 16-51% below HC3 on a 12-point design with one leveraged row, i.e. exactly where robust standard errors are reached for. (robust=True keeps meaning HC0 for compatibility.)

    • HC1: HC0 scaled by n / (n - rank) – the df correction Stata’s robust applies. Fixes the average small-sample bias but ignores WHERE the leverage sits.

    • HC2: per-observation scores scaled by 1 / sqrt(1 - h_i); exactly unbiased under homoskedasticity. Use when leverage is uneven but n is moderate.

    • HC3: scores scaled by 1 / (1 - h_i) (jackknife-like) – the standard recommendation for n <~ 50 or clearly leveraged designs (Long–Ervin 2000); the most conservative of the four.

    All four assume INDEPENDENT observations (none is cluster-robust) and still require the mean/link to be correct (audit G-7). h_i is the observation’s leverage in the converged working weighted regression; under weight_type="frequency" each replicate carries its own leverage, so every estimator equals the expanded data set’s. HC1-3 are refused on a saturated fit (their corrections divide by zero there), and HC2/HC3 on a unit-leverage observation (one fitted by itself alone).

  • theta (float) – the negative-binomial dispersion parameter (family="negativebinomial" only; ignored otherwise). Not estimated from data – pass the value appropriate to your data (e.g. from a prior fit or a method-of-moments estimate); the default 1.0 is a plain placeholder, not a fitted value. WARNING (audit G-4): the MODEL-BASED standard errors assume theta is correct – with the placeholder against true theta = 10, measured SEs were 1.87x the true sampling spread (nominal 95% Wald intervals covering ~100%, tests losing essentially all power). Pass robust=True unless theta is trusted, and if theta came from the same data, remember AIC/BIC do not count it.

Returns:

A GLMResult.

Raises:

PerfectSeparationError – a binomial fit whose fitted probabilities reached exactly 0 or 1, leaving IRLS no finite step. Separation that stalls just short of that returns instead, with GLMResult.separated set, a UserWarning, and z_values / p_values refusing – so treat separated, not the absence of this exception, as the answer to “were my classes separated?” (the exception fires on a continuously coded separation and not on the dummy-coded form of the very same data).

Return type:

GLMResult

class PenalizedResult(coef, intercept, alpha, l1_ratio, n_iter, converged=True, rank=None)[source]

Bases: object

Fitted penalized linear regression.

Parameters:
coef

(p,) coefficients (excluding the intercept).

Type:

numpy.ndarray

intercept

fitted intercept.

Type:

float

alpha

overall penalty strength.

Type:

float

l1_ratio

elastic-net mixing (1 = lasso, 0 = ridge).

Type:

float

n_iter

coordinate-descent iterations (0 for the closed-form ridge).

Type:

int

predict(x)[source]

Predict penalized-regression responses for design rows.

Parameters:

x (ndarray)

Return type:

ndarray

ridge_regression(x, y, alpha=1.0, *, fit_intercept=True)[source]

Ridge (L2-penalised) linear regression in closed form.

Minimises ||y - X b||^2 + alpha ||b||^2; the intercept (if fitted) is not penalised.

Parameters:
Return type:

PenalizedResult

elastic_net(x, y, alpha=1.0, l1_ratio=0.5, *, fit_intercept=True, max_iter=1000, tol=1e-7)[source]

Elastic-net linear regression by cyclic coordinate descent.

Minimises (1/2n) ||y - X b||^2 + alpha ( l1_ratio ||b||_1 + (1 - l1_ratio)/2 ||b||^2 ). l1_ratio = 1 is the lasso (sparse); l1_ratio = 0 is ridge-shaped, but note the 1/(2n) data term means elastic_net(X, y, a, 0.0) equals ridge_regression(X, y, n*a) – the penalties differ by a factor of n when cross-walking alpha (audit G-5).

Parameters:
Return type:

PenalizedResult

lasso(x, y, alpha=1.0, **kw)[source]

Lasso (L1) linear regression – elastic_net() with l1_ratio = 1.

Parameters:
Return type:

PenalizedResult

class RegressionFit(coef, fitted, scale, n_iter, converged=True, rank=None, degenerate_scale=False)[source]

Bases: object

Coefficients + fitted values from robust_regression() / quantile_regression().

Parameters:
degenerate_scale

set by robust_regression() when it detects a sign of its own >50% breakdown point – EITHER the IRLS scale has collapsed (at or within a convergence step of its numerical floor, or far below the raw response’s own spread) with a real share of rows assigned (near) zero weight as a result, OR at least ~44% of the raw response is tied (near enough) to a single value regardless of the fit’s residuals (the shape a capped/censored/zero-inflated measurement produces, which need not collapse the scale at all; see the module comment above _robust_weight_collapse in this file for the full design and, just as importantly, what it still cannot see). Either signal firing means this is not a transient numerical accident, and coef MAY reflect a majority pattern rather than a fully recovered relationship – but may equally reflect a fit that is completely correct, with the flagged pattern arising because the majority genuinely fits (near-)exactly; see robust_regression(). This field, and the UserWarning announcing it at fit time, name the observed CONDITION only, never a verdict: neither can establish whether a real minority signal was discarded or a genuine minority of contamination was correctly rejected (both look identical from this pattern alone), and – a further, structural limitation, not merely an unresolved ambiguity – a majority and minority that follow two genuinely different relationships (as opposed to sharing one tied value) is largely invisible to this mechanism in the first place unless the majority’s own fit is exceptionally tight; see the module comment for exactly how far that reaches and where it stops. Always False for quantile_regression(), which has no comparable mechanism.

Type:

bool

predict(x)[source]

Predict fitted regression values for design rows.

Parameters:

x (ndarray)

Return type:

ndarray

robust_regression(x, y, *, method='huber', c=None, max_iter=100, tol=1e-8)[source]

Robust (M-estimator) linear regression by IRLS with a robust scale.

Down-weights observations with large residuals so a few outliers cannot dominate the fit. huber uses the Huber weight (tuning c = 1.345 for 95% Gaussian efficiency); tukey uses the redescending Tukey biweight (c = 4.685), which rejects gross outliers entirely.

The scale estimate (median absolute residual deviation / 0.6745) is only ever as robust as any MAD-based estimator is: past its own breakdown point – half or more of the rows sharing (near enough) the same residual under the current fit, an ordinary shape for zero-inflated or capped/rounded data and not a contrived one – it is exactly zero by construction and floors at 1e-8 (or, as of audit R-3, is treated the same way once it is merely far smaller than the response’s own spread, even without literally reaching that floor). That is an inherent limit of this estimator family, not something iterating further or switching method fixes: tukey’s redescending weight then hard-zeros essentially every row (raised below), and huber’s never-quite-zero weight instead settles – sometimes only after many iterations – on a near-zero-coefficient fit matching the majority tie that reports converged=True. A capped or censored measurement breaks down the same way without ever collapsing the scale estimate itself: a shared ceiling ties rows’ raw responses together but not their residuals once the fit carries any slope, so this is instead read directly off the response (see the module comment above _robust_weight_collapse in this file for the full two-signal detection design, including what audit R-3 closed and what it plainly did not: most notably, a majority and a minority that follow two genuinely DIFFERENT relationships – as opposed to sharing one tied value – is only ever PARTIALLY visible to this mechanism, and a location-only split with no covariate to anchor the response’s own spread is not visible to it at all). All of these are disclosed rather than left looking like an ordinary converged fit: see RegressionFit.degenerate_scale and the UserWarning issued at fit time. What is NOT, and cannot be, disclosed: WHICH of “a real minority signal was discarded,” “a minority of contamination was correctly rejected,” or “nothing was discarded at all, because the majority is genuinely this exact” produced a flagged pattern – all three leave the same fingerprint, so the warning names the observed condition, never a verdict on which case it is.

A sharper, method='huber'-specific gap (GAP C in the module comment, campaign nine, D-0209): huber’s weight never fully crushes, so a real minority following a different relationship can be systematically down-weighted enough to bias coef toward the majority WITHOUT tripping either disclosure signal at all – confirmed common (up to 100% of seeds at a 35% minority fraction, at everyday noise levels), not a rare edge case. Raising max_iter past a RuntimeError("... failed to converge ...") is what usually exposes it: that default failure is a real but incomplete guard, not a guarantee, so a non-convergence error under method='huber' is worth treating as a possible sign of exactly this shape rather than simply raising max_iter and moving on. Switching to method='tukey' is not a reliable workaround for this specific shape either outside a no-intercept design.

Raises:

RuntimeError – every observation lands at zero weight (the message names the breakdown pattern above when one was detected), coefficients become non-finite, or IRLS does not converge in max_iter iterations.

Parameters:
Return type:

RegressionFit

quantile_regression(x, y, tau=0.5, *, max_iter=200, tol=1e-7, eps=1e-6)[source]

Linear quantile regression: the conditional tau-quantile by IRLS on the check loss.

Minimises the pinball loss sum rho_tau(y - X b) via iteratively reweighted least squares with weights tau / |r| for positive residuals and (1 - tau) / |r| for negative ones (a smoothed Newton scheme; eps floors |r| for stability).

Parameters:
Return type:

RegressionFit