# Family


Abstract family defining the response distribution and link function.


Usage

``` python
Family()
```


A [Family](Family.md#whittaker.Family) encapsulates everything the P-IRLS fitting loop needs to know about the conditional distribution of the response `y` given the linear predictor `eta`. Every GLM and GAM family in Whittaker belongs to the exponential dispersion family, and is fully characterized by three ingredients:

1.  The **link function** `g`, relating the mean `mu` to the linear predictor `eta`: `eta = g(mu)`, along with its inverse and derivative.
2.  The **variance function** `V(mu)`, relating the variance of the response to its mean: `Var(Y) = phi * V(mu)`, where `phi` is the dispersion (scale) parameter.
3.  The **deviance** and **log-likelihood**, which quantify goodness of fit and are used by [GAM.fit()](GAM.md#whittaker.GAM.fit) for smoothing parameter selection (GCV/REML) and by [GAM.summary()](GAM.md#whittaker.GAM.summary) for reporting.

Subclasses must implement all abstract methods (`link`, `link_inverse`, `link_derivative`, `variance`, `deviance`, `log_likelihood`, `simulate`). The link function and its inverse/derivative are used by the P-IRLS algorithm to form pseudo-data (the working response `z`) and working weights `W` at each iteration. Families whose loss does not fit the standard GLM deviance framework (e.g. [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily), [CoxPH](CoxPH.md#whittaker.CoxPH), [OrderedCategorical](OrderedCategorical.md#whittaker.OrderedCategorical), [Multinomial](Multinomial.md#whittaker.Multinomial)) may instead override `irls_update` to supply `z` and `W` directly.

Whittaker ships with the following concrete families:

- [Gaussian](Gaussian.md#whittaker.Gaussian) -- identity link, constant variance; the default family for continuous, unbounded responses.
- [Poisson](Poisson.md#whittaker.Poisson) -- log link, `V(mu) = mu`; for count data.
- [Binomial](Binomial.md#whittaker.Binomial) -- logit link, `V(mu) = mu(1-mu)`; for binary or proportion responses.
- [Gamma](Gamma.md#whittaker.Gamma) -- log link, `V(mu) = mu^2`; for positive, right-skewed continuous data.
- [NegativeBinomial](NegativeBinomial.md#whittaker.NegativeBinomial) -- log link, `V(mu) = mu + mu^2/theta`; for overdispersed counts.
- [Beta](Beta.md#whittaker.Beta) -- logit link; for proportions strictly between 0 and 1.
- [Tweedie](Tweedie.md#whittaker.Tweedie) / [TweedieEstimated](TweedieEstimated.md#whittaker.TweedieEstimated) (via [tw()](tw.md#whittaker.tw)) -- log link, `V(mu) = mu^p`; for compound Poisson-Gamma data with a point mass at zero (e.g. insurance claims).
- [InverseGaussian](InverseGaussian.md#whittaker.InverseGaussian) -- log link, `V(mu) = mu^3`; for positive, heavy-tailed continuous data.
- [CoxPH](CoxPH.md#whittaker.CoxPH) -- proportional hazards partial likelihood; for survival/time-to-event data.
- [OrderedCategorical](OrderedCategorical.md#whittaker.OrderedCategorical) -- cumulative logit (proportional odds); for ordinal responses.
- [Multinomial](Multinomial.md#whittaker.Multinomial) -- baseline-category logit; for unordered categorical responses.
- [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily) -- Extended Log-F (ELF) smooth pinball loss; for quantile regression.

For distributional regression, where more than one parameter of the response distribution (e.g. both the mean and the scale) is modeled by its own smooth predictor, see [GAMLSSFamily](GAMLSSFamily.md#whittaker.GAMLSSFamily) and its concrete subclasses ([GaussianLS](GaussianLS.md#whittaker.GaussianLS), [GammaLS](GammaLS.md#whittaker.GammaLS), [BetaLS](BetaLS.md#whittaker.BetaLS), [ZeroInflatedPoisson](ZeroInflatedPoisson.md#whittaker.ZeroInflatedPoisson), [ZeroInflatedNegativeBinomial](ZeroInflatedNegativeBinomial.md#whittaker.ZeroInflatedNegativeBinomial)).


## Examples

Families are passed to [GAM](GAM.md#whittaker.GAM) via the `family` argument; they are rarely instantiated by users beyond that.


``` python
import numpy as np
import whittaker as wk

rng = np.random.default_rng(0)
n = 200
x = np.linspace(0, 2 * np.pi, n)
y = rng.poisson(np.exp(0.5 * np.sin(x)))

data = {"x": x, "y": y}

# Any concrete Family subclass can be passed to GAM
model = wk.GAM("y ~ s(x)", family=wk.Poisson())
model.fit(data, method="REML")
print(model.summary())
```


    GAM fit summary
    ============================================================
    Formula:    y ~ s(x)
    Family:     Poisson(link='log')
    Inference:  REML
    Observations: 200
    Coefficients: 10

    Parametric coefficients:
      Term                       Estimate    Std.Err    z value    p-value
      ------------------------ ---------- ---------- ---------- ----------
      (Intercept)                  0.2113     0.0657      3.214    0.00131

    Approximate significance of smooth terms:
      Term                        EDF Ref.df     Chi.sq    p-value
      ------------------------ ------ ------ ---------- ----------
      s(x)                       4.18      5     37.260  5.312e-07

    Total EDF:  5.18
    Scale est:  1.000000
    Deviance:   187.0851
    Null dev:   230.5815
    Dev. expl:  18.9%
    GCV score:  0.985797
    AIC:        554.30
    BIC:        571.37


## Attributes

| Name | Description |
|----|----|
| [scale_known](#scale_known) | Whether the dispersion (scale) parameter is fixed rather than estimated. |

------------------------------------------------------------------------


### scale_known


Whether the dispersion (scale) parameter is fixed rather than estimated.


`scale_known: bool`


Returns `True` for families with no free dispersion parameter (e.g. [Binomial](Binomial.md#whittaker.Binomial), [Poisson](Poisson.md#whittaker.Poisson)), in which case the scale is always `1` and is not estimated during fitting. Returns `False` (the default) for families such as [Gaussian](Gaussian.md#whittaker.Gaussian) and [Gamma](Gamma.md#whittaker.Gamma), whose scale is estimated from the data.


## Methods

| Name | Description |
|----|----|
| [deviance()](#deviance) | Total (unscaled) deviance. |
| [initialize()](#initialize) | Compute starting values for `mu` given the observed response `y`. |
| [irls_update()](#irls_update) | Custom IRLS pseudo-response and working weights. |
| [link()](#link) | Apply the link function g(\mu), mapping the conditional mean to the linear predictor. |
| [link_derivative()](#link_derivative) | Derivative of the link function, g'(\mu) = d\eta/d\mu. |
| [link_inverse()](#link_inverse) | Apply the inverse link g^{-1}(\eta), mapping the linear predictor back to the mean. |
| [log_lik_pointwise()](#log_lik_pointwise) | Per-observation log-likelihood contributions \ell_i(y_i; \mu_i, \phi). |
| [log_likelihood()](#log_likelihood) | Log-likelihood \ell(y; \mu, \phi) evaluated at the given scale parameter \phi. |
| [simulate()](#simulate) | Simulate response values from the distribution. |
| [unit_deviance()](#unit_deviance) | Per-observation deviance contributions d_i, before summing over observations. |
| [variance()](#variance) | Variance function V(\mu) relating the response variance to its mean. |

------------------------------------------------------------------------


### deviance()


Total (unscaled) deviance.


Usage

``` python
deviance(
    y,
    mu,
    *,
    weights=None,
)
```


D(y, \hat\mu) = 2 \sum_i \[\ell(y_i; y_i) - \ell(y_i; \hat\mu_i)\].

The deviance measures the discrepancy between the fitted model and a saturated model that fits the data exactly. It is used by [GAM.fit()](GAM.md#whittaker.GAM.fit) for smoothing parameter selection (GCV) and by [GAM.summary()](GAM.md#whittaker.GAM.summary) for goodness-of-fit reporting.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.

`weights: NDArray | None = None`  
Optional prior weights, shape `(n,)`. When given, the deviance is \sum_i w_i d_i rather than \sum_i d_i.


#### Returns


`float`  
The total deviance.


------------------------------------------------------------------------


### initialize()


Compute starting values for `mu` given the observed response `y`.


Usage

``` python
initialize(y)
```


Called once before the first P-IRLS iteration to seed the mean. The default implementation returns `y` unchanged, which is appropriate for the [Gaussian](Gaussian.md#whittaker.Gaussian) family with its identity link. Families with non-identity links or constrained means (e.g. [Poisson](Poisson.md#whittaker.Poisson), [Binomial](Binomial.md#whittaker.Binomial), [Gamma](Gamma.md#whittaker.Gamma)) override this to nudge `y` into the valid range and avoid numerical issues (such as `log(0)`) on the first iteration.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.


#### Returns


`NDArray`  
Starting values for `mu`, shape `(n,)`.


------------------------------------------------------------------------


### irls_update()


Custom IRLS pseudo-response and working weights.


Usage

``` python
irls_update(
    y,
    mu,
    eta,
)
```


Override in families whose loss does not fit the standard GLM deviance framework (e.g. [QuantileFamily](QuantileFamily.md#whittaker.QuantileFamily), [CoxPH](CoxPH.md#whittaker.CoxPH), [OrderedCategorical](OrderedCategorical.md#whittaker.OrderedCategorical), [Multinomial](Multinomial.md#whittaker.Multinomial)), to supply the working response `z` and working weights `W` directly rather than deriving them from `link`, `link_derivative`, and `variance`. Returning `None` (the default) tells the P-IRLS loop to fall back to the standard GLM formula `z = eta + (y - mu) * link_derivative(mu)`, `W = 1 / (link_derivative(mu)^2 * variance(mu))`.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`mu: NDArray`  
Current fitted mean values, shape `(n,)`.

`eta: NDArray`  
Current linear predictor values, shape `(n,)`.


#### Returns


`tuple of NDArray, or None`  
`(z, W)`, the pseudo-response and diagonal working-weight vector, each shape `(n,)`, or `None` to use the default GLM formula.


------------------------------------------------------------------------


### link()


Apply the link function g(\mu), mapping the conditional mean to the linear predictor.


Usage

``` python
link(mu)
```


The link function defines the relationship \eta = g(\mu) between the mean of the response and the linear predictor. Each family provides a canonical or default link; for example, the identity link for Gaussian, the log link for Poisson, and the logit link for Binomial.


#### Parameters


`mu: NDArray`  
Conditional mean values \mu, shape `(n,)`. Must lie in the valid range for the family (e.g., \mu \> 0 for Poisson, 0 \< \mu \< 1 for Binomial).


#### Returns


`NDArray`  
Linear predictor values \eta = g(\mu), shape `(n,)`.


------------------------------------------------------------------------


### link_derivative()


Derivative of the link function, g'(\mu) = d\eta/d\mu.


Usage

``` python
link_derivative(mu)
```


Used by the P-IRLS fitting loop to form the working response `z` and working weights `W` at each iteration, via a first-order (delta-method) linearization of the link function around the current fit.


#### Parameters


`mu: NDArray`  
Conditional mean values \mu, shape `(n,)`.


#### Returns


`NDArray`  
Derivative values g'(\mu), shape `(n,)`.


------------------------------------------------------------------------


### link_inverse()


Apply the inverse link g^{-1}(\eta), mapping the linear predictor back to the mean.


Usage

``` python
link_inverse(eta)
```


This is the transformation applied to the fitted linear predictor to recover fitted values `mu` on the response scale, e.g. after [GAM.predict()](GAM.md#whittaker.GAM.predict) computes `eta`.


#### Parameters


`eta: NDArray`  
Linear predictor values \eta, shape `(n,)`.


#### Returns


`NDArray`  
Conditional mean values \mu = g^{-1}(\eta), shape `(n,)`.


------------------------------------------------------------------------


### log_lik_pointwise()


Per-observation log-likelihood contributions \ell_i(y_i; \mu_i, \phi).


Usage

``` python
log_lik_pointwise(
    y,
    mu,
    scale,
    *,
    weights=None,
)
```


Returns the same values that `log_likelihood` sums, without summing, as a 1-D array. Used by [GAM.loo()](GAM.md#whittaker.GAM.loo) to compute PSIS-LOO cross-validation.

The default implementation calls `log_likelihood` on each observation individually. Key families override this with a vectorized implementation.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.

`scale: float`  
Dispersion (scale) parameter \phi.

`weights: NDArray | None = None`  
Optional prior weights, shape `(n,)`.


#### Returns


`NDArray`  
Per-observation log-likelihood values, shape `(n,)`.


------------------------------------------------------------------------


### log_likelihood()


Log-likelihood \ell(y; \mu, \phi) evaluated at the given scale parameter \phi.


Usage

``` python
log_likelihood(
    y,
    mu,
    scale,
    *,
    weights=None,
)
```


Used by [GAM.fit()](GAM.md#whittaker.GAM.fit) for REML/ML-based smoothing parameter selection and by [GAM.summary()](GAM.md#whittaker.GAM.summary) for reporting AIC and related fit statistics.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.

`scale: float`  
Dispersion (scale) parameter \phi.

`weights: NDArray | None = None`  
Optional prior weights, shape `(n,)`. When given, the log-likelihood is \sum_i w_i \ell_i rather than \sum_i \ell_i.


#### Returns


`float`  
The total log-likelihood.


------------------------------------------------------------------------


### simulate()


Simulate response values from the distribution.


Usage

``` python
simulate(
    mu,
    scale,
    rng,
)
```


#### Parameters


`mu: NDArray`  
Mean (fitted values), shape `(n,)`.

`scale: float`  
Estimated scale parameter φ.

`rng: np.random.Generator`  
A `numpy.random.Generator` instance.


#### Returns


`NDArray`  
Simulated response values, shape `(n,)`.


------------------------------------------------------------------------


### unit_deviance()


Per-observation deviance contributions d_i, before summing over observations.


Usage

``` python
unit_deviance(
    y,
    mu,
)
```


The total deviance returned by `deviance` is \sum_i d_i (or the weighted sum). The default implementation here is `(y - mu)^2`, appropriate for the Gaussian family; concrete families with a different deviance formula override this method.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`mu: NDArray`  
Fitted conditional mean values, shape `(n,)`.


#### Returns


`NDArray`  
Per-observation deviance contributions d_i, shape `(n,)`.


------------------------------------------------------------------------


### variance()


Variance function V(\mu) relating the response variance to its mean.


Usage

``` python
variance(mu)
```


The conditional variance of the response is \operatorname{Var}(Y) = \phi \\ V(\mu), where \phi is the dispersion (scale) parameter. For the Gaussian family this is constant (`1`); for Poisson it is `mu`; for Gamma it is `mu^2`, etc. `variance` is used by P-IRLS to form the working weights.


#### Parameters


`mu: NDArray`  
Conditional mean values \mu, shape `(n,)`.


#### Returns


`NDArray`  
Variance function values V(\mu), shape `(n,)`.