GLM Families API

The GLM module provides distribution family abstractions for P-spline smoothing beyond Gaussian responses, using Iteratively Reweighted Least Squares (IRLS).

Family Protocol

Bases: Protocol

              flowchart TD
              psplines.glm.Family[Family]

              

              click psplines.glm.Family href "" "psplines.glm.Family"
            

Protocol for GLM family implementations.

is_gaussian property

is_gaussian: bool

Whether this is the Gaussian (identity link) family.

initialize

initialize(
    y: NDArray, **kwargs: object
) -> tuple[NDArray, NDArray, NDArray]

Return (eta, mu, W_diag) starting values.

Source code in src/psplines/glm.py

def initialize(
    self, y: NDArray, **kwargs: object
) -> tuple[NDArray, NDArray, NDArray]:
    """Return (eta, mu, W_diag) starting values."""
    ...

working_response

working_response(
    y: NDArray, eta: NDArray, mu: NDArray
) -> NDArray

Compute working dependent variable z = eta + W^{-1}(y - mu).

Source code in src/psplines/glm.py

def working_response(self, y: NDArray, eta: NDArray, mu: NDArray) -> NDArray:
    """Compute working dependent variable z = eta + W^{-1}(y - mu)."""
    ...

working_weights

working_weights(mu: NDArray, **kwargs: object) -> NDArray

Compute diagonal of IRLS weight matrix W.

Source code in src/psplines/glm.py

def working_weights(self, mu: NDArray, **kwargs: object) -> NDArray:
    """Compute diagonal of IRLS weight matrix W."""
    ...

deviance

deviance(y: NDArray, mu: NDArray) -> float

Compute model deviance.

Source code in src/psplines/glm.py

def deviance(self, y: NDArray, mu: NDArray) -> float:
    """Compute model deviance."""
    ...

inverse_link

inverse_link(eta: NDArray) -> NDArray

Map linear predictor to mean: mu = h(eta).

Source code in src/psplines/glm.py

def inverse_link(self, eta: NDArray) -> NDArray:
    """Map linear predictor to mean: mu = h(eta)."""
    ...

phi

phi(deviance: float, n: int, ed: float) -> float

Compute scale parameter phi.

Source code in src/psplines/glm.py

def phi(self, deviance: float, n: int, ed: float) -> float:
    """Compute scale parameter phi."""
    ...

GaussianFamily

Gaussian family with identity link.

This wraps the existing Gaussian P-spline behaviour so the IRLS loop converges in a single iteration.

PoissonFamily

Poisson family with log link (§2.12.1, eq. 2.18–2.22).

Parameters:

• offset (NDArray or None, default: None ) –

  Log-scale exposure offset. When provided, mu = exp(eta) * exp(offset), i.e. eta_total = B @ alpha + offset.

Source code in src/psplines/glm.py

def __init__(self, offset: NDArray | None = None) -> None:
    self.offset = offset

BinomialFamily

Binomial family with logit link (§2.12.2, eq. 2.23–2.25).

Parameters:

• trials (NDArray or None, default: None ) –

  Number of trials per observation. Defaults to 1 (Bernoulli).

Source code in src/psplines/glm.py

def __init__(self, trials: NDArray | None = None) -> None:
    self.trials = trials

get_family

Resolve a family string or instance into a Family object.

Parameters:

• family (str or Family) –

  One of "gaussian", "poisson", "binomial", or a Family instance.

• trials (NDArray, default: None ) –

  Binomial trials vector.

• offset (NDArray, default: None ) –

  Poisson offset (log-exposure).

Returns:

• Family –

  Resolved family instance.

Source code in src/psplines/glm.py

def get_family(
    family: str | Family,
    trials: NDArray | None = None,
    offset: NDArray | None = None,
) -> Family:
    """
    Resolve a family string or instance into a Family object.

    Parameters
    ----------
    family : str or Family
        One of "gaussian", "poisson", "binomial", or a Family instance.
    trials : NDArray, optional
        Binomial trials vector.
    offset : NDArray, optional
        Poisson offset (log-exposure).

    Returns
    -------
    Family
        Resolved family instance.
    """
    if isinstance(family, str):
        name = family.lower()
        if name == "gaussian":
            return GaussianFamily()
        if name == "poisson":
            return PoissonFamily(offset=offset)
        if name == "binomial":
            return BinomialFamily(trials=trials)
        raise ValueError(
            f"Unknown family '{family}'. Choose from: 'gaussian', 'poisson', 'binomial'"
        )
    return family

Overview

The GLM framework extends P-splines to non-Gaussian responses via IRLS. Instead of solving the normal equations directly, the fitting loop iterates:

1. Compute working response \(z = \eta + W^{-1}(y - \mu)\)
2. Compute working weights \(W\) from the current \(\mu\)
3. Solve \((B'WB + \lambda D'D)\alpha = B'Wz\)
4. Update \(\eta = B\alpha\), \(\mu = h(\eta)\)
5. Repeat until convergence

Each family defines the link function \(h\), the working weights, and the deviance.

Supported Families

Family	Link	\(\mu = h(\eta)\)	Working weights \(W\)	Scale \(\phi\)
Gaussian	Identity	\(\eta\)	\(1\)	Estimated
Poisson	Log	\(\exp(\eta)\)	\(\text{diag}(\mu)\)	1
Binomial	Logit	\(t \cdot \text{sigmoid}(\eta)\)	\(\text{diag}(\mu(1-\pi))\)	1

Usage with PSpline

Pass a family string (or instance) to the PSpline constructor:

from psplines import PSpline

# Poisson (count data)
spline = PSpline(x, counts, family="poisson")

# Poisson with exposure offset (rate modeling)
spline = PSpline(x, counts, family="poisson", offset=np.log(exposure))

# Binomial — Bernoulli (binary response)
spline = PSpline(x, binary_y, family="binomial")

# Binomial — grouped (y successes out of t trials)
spline = PSpline(x, successes, family="binomial", trials=trials_vec)

Custom Families

Any object satisfying the Family protocol can be passed directly:

from psplines.glm import Family

class MyFamily:
    """Custom family implementing the Family protocol."""

    @property
    def is_gaussian(self) -> bool:
        return False

    def initialize(self, y, **kwargs):
        ...

    def working_response(self, y, eta, mu):
        ...

    def working_weights(self, mu, **kwargs):
        ...

    def deviance(self, y, mu):
        ...

    def inverse_link(self, eta):
        ...

    def phi(self, deviance, n, ed):
        ...

spline = PSpline(x, y, family=MyFamily())