# Response families
A GAM's **family** specifies two things: the conditional distribution of the response y given the predictors, and the **link function** g that relates the conditional mean \mu to the linear predictor \eta:
g(\mu) = \eta = \beta_0 + f_1(x_1) + f_2(x_2) + \cdots
The choice of family determines how predictions are mapped back to the response scale, how variance changes with the mean, and which loss function (deviance) is optimized during fitting. Getting the family right is essential: a misspecified family can bias estimates, distort confidence intervals, and produce nonsensical predictions.
# The exponential family framework
All families in Whittaker belong to the exponential dispersion family, characterized by a density of the form:
f(y \mid \theta, \phi) = a(y, \phi)\\\exp\\\left\[\frac{y\theta - b(\theta)}{\phi}\right\]
where \theta is the **canonical parameter**, \phi is the **dispersion parameter**, and b(\cdot) is the **cumulant function**. From this formulation, four key ingredients follow:
- **Mean**: \mu = b'(\theta).
- **Variance function**: \text{Var}(Y) = \phi\\ V(\mu), where V(\mu) = b''(\theta). The shape of V(\mu) is what distinguishes one family from another.
- **Link function**: g(\mu) = \eta. The **canonical link** sets g(\mu) = \theta (alternative links are also supported).
- **Deviance**: D(y, \hat\mu) = 2\phi \sum_i \bigl\[\ell(y_i; y_i) - \ell(y_i; \hat\mu_i)\bigr\], where \ell is the log-likelihood. Deviance is the quantity minimized during P-IRLS fitting (see [Model fitting](fitting.md)).
> **Tip: Canonical vs. non-canonical links**
>
> Using the canonical link simplifies the score equations and guarantees concavity of the log-likelihood. Non-canonical links (e.g., `log` for Gamma instead of `inverse`) are often preferred for interpretability but may require more P-IRLS iterations to converge.
# Quick reference table
The table below summarizes every family available in Whittaker. GAMLSS families (for distributional regression) are covered in a [separate page](gamlss.md).
| Family | String | Default link | V(\mu) | Typical use |
|----|----|----|----|----|
| [Gaussian()](../reference/Gaussian.md#whittaker.Gaussian) | `"gaussian"` | identity | 1 | Continuous, unbounded |
| [Poisson()](../reference/Poisson.md#whittaker.Poisson) | `"poisson"` | log | \mu | Event counts |
| [Binomial()](../reference/Binomial.md#whittaker.Binomial) | `"binomial"` | logit | \mu(1-\mu) | Binary / proportions |
| [Gamma()](../reference/Gamma.md#whittaker.Gamma) | `"gamma"` | inverse | \mu^2 | Positive continuous (right-skewed) |
| [NegativeBinomial()](../reference/NegativeBinomial.md#whittaker.NegativeBinomial) | `"nb"` | log | \mu + \mu^2/\theta | Overdispersed counts |
| [Beta()](../reference/Beta.md#whittaker.Beta) | `"beta"` | logit | \mu(1-\mu)/(1+\phi) | Proportions on (0,1) |
| [Tweedie()](../reference/Tweedie.md#whittaker.Tweedie) | `"tweedie"` | log | \mu^p | Mixed zero/continuous |
| [InverseGaussian()](../reference/InverseGaussian.md#whittaker.InverseGaussian) | `"inverse.gaussian"` | 1/\mu^2 | \mu^3 | Positive, heavy-tailed |
| [CoxPH()](../reference/CoxPH.md#whittaker.CoxPH) | `"coxph"` | log | | Survival / hazard |
# Specifying a family
There are two ways to set the family: a string shorthand or a family object. Use the string when the default link is what you want. Use the object when you need a non-default link or want to set initial parameter values.
``` python
# String shorthand
model = wk.GAM("y ~ s(x)", family="gaussian")
# Family object
model = wk.GAM("y ~ s(x)", family=wk.Gaussian())
```
# Gaussian: continuous responses
The Gaussian family is the default. It assumes the response is continuous and unbounded, with constant variance:
- **Link**: g(\mu) = \mu (identity).
- **Variance function**: V(\mu) = 1, so \text{Var}(Y) = \phi.
- **Deviance**: D = \sum_i (y_i - \hat\mu_i)^2 (residual sum of squares).
This is the classical regression setting: the model minimizes penalized least squares and converges in a single step of P-IRLS.
**Supported links**: `identity` (default), `log`, `inverse`.
## Example: noisy sine curve
``` python
import numpy as np
import whittaker as wk
import altair as alt
rng = np.random.default_rng(0)
n = 200
x = np.linspace(0, 2 * np.pi, n)
y = np.sin(x) + rng.normal(0, 0.3, n)
data = {"x": x, "y": y}
model_gauss = wk.GAM("y ~ s(x)", family=wk.Gaussian())
model_gauss.fit(data, method="REML")
print(model_gauss.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Gaussian(link='identity')
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.0046 0.0205 0.224 0.8233
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 7.31 8 1279.244 < 1e-16
Total EDF: 8.31
Deviance: 16.0702
Null dev: 123.9418
Dev. expl: 87.0%
GCV score: 0.087470
Scale est: 0.083835
AIC: 80.11
BIC: 107.52
Because the identity link maps the linear predictor directly to the mean, predictions are on the original scale without any back-transformation.
## Visualizing the Gaussian fit
``` python
x_grid = np.linspace(0, 2 * np.pi, 300)
preds = model_gauss.predict({"x": x_grid}, se=True)
z = 1.96
lower = preds.values - z * preds.se
upper = preds.values + z * preds.se
points = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y[i])} for i in range(n)]}
).mark_circle(size=15, opacity=0.3, color="steelblue").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="y"),
)
fit_data = [
{"x": float(x_grid[i]), "fit": float(preds.values[i]),
"lower": float(lower[i]), "upper": float(upper[i])}
for i in range(len(x_grid))
]
line = alt.Chart({"values": fit_data}).mark_line(
color="firebrick", strokeWidth=2
).encode(x="x:Q", y="fit:Q")
band = alt.Chart({"values": fit_data}).mark_area(
opacity=0.2, color="firebrick"
).encode(x="x:Q", y="lower:Q", y2="upper:Q")
(band + points + line).properties(
width="container", height=300,
title="Gaussian GAM: y ~ s(x)"
)
```
# Poisson: count responses
The Poisson family is the standard choice for non-negative integer counts. Its key assumption is that the variance equals the mean (equidispersion).
- **Link**: g(\mu) = \log(\mu) (log link).
- **Variance function**: V(\mu) = \mu.
- **Deviance**: D = 2\sum_i \bigl\[y_i \log(y_i / \hat\mu_i) - (y_i - \hat\mu_i)\bigr\].
The log link ensures that predicted counts are always non-negative.
**Supported links**: `log` (default), `identity`, `sqrt`.
> **Tip: Exposure and offsets**
>
> When counts are observed over varying durations or areas, include an offset term to model the **rate** rather than the raw count:
>
> ``` python
> model = wk.GAM("count ~ s(x) + offset(log_exposure)", family="poisson")
> ```
>
> Here `log_exposure` is the natural log of the exposure variable.
## Example: species counts along a transect
``` python
rng = np.random.default_rng(1)
n = 250
x = np.linspace(0, 6, n)
eta = 0.8 + 1.2 * np.sin(x) # log-scale linear predictor
mu = np.exp(eta) # true rate
y_counts = rng.poisson(mu).astype(float)
data_pois = {"x": x, "y": y_counts}
model_pois = wk.GAM("y ~ s(x, k=15)", family=wk.Poisson())
model_pois.fit(data_pois, method="REML")
print(model_pois.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x, k=15)
Family: Poisson(link='log')
Observations: 250
Coefficients: 15
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.8603 0.0479 17.973 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x, k=15) 6.63 7 338.219 < 1e-16
Total EDF: 7.63
Deviance: 263.2131
Null dev: 704.1532
Dev. expl: 62.6%
GCV score: 1.120213
Scale est: 1.000000
AIC: 897.95
BIC: 924.83
## Visualizing the Poisson fit
``` python
x_grid = np.linspace(0, 6, 300)
preds_pois = model_pois.predict({"x": x_grid}, se=True)
# Confidence intervals on the response (count) scale
z = 1.96
lower_pois = np.exp(np.log(preds_pois.values) - z * preds_pois.se)
upper_pois = np.exp(np.log(preds_pois.values) + z * preds_pois.se)
points_pois = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_counts[i])} for i in range(n)]}
).mark_circle(size=15, opacity=0.3, color="teal").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="Count"),
)
fit_data_pois = [
{"x": float(x_grid[i]), "fit": float(preds_pois.values[i]),
"lower": float(lower_pois[i]), "upper": float(upper_pois[i])}
for i in range(len(x_grid))
]
line_pois = alt.Chart({"values": fit_data_pois}).mark_line(
color="darkorange", strokeWidth=2
).encode(x="x:Q", y="fit:Q")
band_pois = alt.Chart({"values": fit_data_pois}).mark_area(
opacity=0.15, color="darkorange"
).encode(x="x:Q", y="lower:Q", y2="upper:Q")
true_pois = [{"x": float(x_grid[i]), "mu": float(np.exp(0.8 + 1.2 * np.sin(x_grid[i])))}
for i in range(len(x_grid))]
true_line_pois = alt.Chart({"values": true_pois}).mark_line(
color="gray", strokeDash=[4, 4], strokeWidth=1.5
).encode(x="x:Q", y="mu:Q")
(band_pois + points_pois + line_pois + true_line_pois).properties(
width="container", height=300,
title="Poisson GAM: count ~ s(x)"
)
```
The orange curve is the estimated rate \hat\mu(x), the gray dashed line is the true generating rate, and each blue-green point is an observed count.
# Binomial: binary and proportion responses
The Binomial family handles two types of response: binary outcomes (0/1) and proportions with a known denominator. The logit link maps probabilities from (0,1) to the real line.
- **Link**: g(\mu) = \log\\\bigl(\mu / (1 - \mu)\bigr) (logit).
- **Variance function**: V(\mu) = \mu(1 - \mu).
- **Deviance**: D = -2\sum_i \bigl\[y_i \log(\hat\mu_i) + (1 - y_i)\log(1 - \hat\mu_i)\bigr\] (for binary data).
**Supported links**: `logit` (default), `probit`, `cloglog`, `cauchit`.
## Example: probability of occurrence
``` python
rng = np.random.default_rng(2)
n = 300
x = np.linspace(-3, 3, n)
eta = -1.0 + 1.5 * x - 0.3 * x**2 # logit-scale predictor
prob = 1 / (1 + np.exp(-eta)) # true probability
y_bin = rng.binomial(1, prob).astype(float)
data_bin = {"x": x, "y": y_bin}
model_bin = wk.GAM("y ~ s(x)", family=wk.Binomial())
model_bin.fit(data_bin, method="REML")
print(model_bin.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Binomial(link='logit')
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) -1.5512 0.3410 -4.549 5.39e-06
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 2.53 3 61.483 2.834e-13
Total EDF: 3.53
Deviance: 231.7358
Null dev: 390.8857
Dev. expl: 40.7%
GCV score: 0.790951
Scale est: 1.000000
AIC: 238.79
BIC: 251.86
## Visualizing the Binomial fit
``` python
x_grid = np.linspace(-3, 3, 300)
preds_bin = model_bin.predict({"x": x_grid}, se=True)
# Confidence intervals via the logit transform
logit_hat = np.log(preds_bin.values / (1 - preds_bin.values))
lower_bin = 1 / (1 + np.exp(-(logit_hat - 1.96 * preds_bin.se)))
upper_bin = 1 / (1 + np.exp(-(logit_hat + 1.96 * preds_bin.se)))
points_bin = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_bin[i])} for i in range(n)]}
).mark_circle(size=12, opacity=0.15, color="purple").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="P(y = 1)"),
)
fit_data_bin = [
{"x": float(x_grid[i]), "fit": float(preds_bin.values[i]),
"lower": float(lower_bin[i]), "upper": float(upper_bin[i])}
for i in range(len(x_grid))
]
line_bin = alt.Chart({"values": fit_data_bin}).mark_line(
color="darkviolet", strokeWidth=2
).encode(x="x:Q", y="fit:Q")
band_bin = alt.Chart({"values": fit_data_bin}).mark_area(
opacity=0.15, color="darkviolet"
).encode(x="x:Q", y="lower:Q", y2="upper:Q")
(band_bin + points_bin + line_bin).properties(
width="container", height=300,
title="Binomial GAM: binary outcome ~ s(x)"
)
```
> **Note: Proportion data**
>
> For grouped binomial data (successes out of n trials), use the `cbind` syntax:
>
> ``` python
> model = wk.GAM("cbind(successes, failures) ~ s(x)", family="binomial")
> ```
>
> This is more efficient than expanding to individual binary rows and correctly accounts for the denominator.
# Gamma: positive continuous responses
The Gamma family is appropriate for strictly positive, right-skewed continuous data where the variance increases with the mean. Insurance claims, precipitation amounts, and response times are classic use cases.
- **Link**: g(\mu) = 1/\mu (inverse, canonical) or g(\mu) = \log(\mu) (log).
- **Variance function**: V(\mu) = \mu^2.
- **Deviance**: D = 2\sum_i \bigl\[-\log(y_i / \hat\mu_i) + (y_i - \hat\mu_i)/\hat\mu_i\bigr\].
The Gamma family uses a **log link** by default, which guarantees positive fitted values.
## Example: reaction time data
``` python
rng = np.random.default_rng(3)
n = 250
x = np.linspace(0.5, 5, n)
mu = np.exp(1.0 + 0.4 * np.sin(2 * x)) # true mean (always positive)
shape = 10.0
y_gamma = rng.gamma(shape, scale=mu / shape, size=n)
data_gamma = {"x": x, "y": y_gamma}
model_gamma = wk.GAM("y ~ s(x)", family=wk.Gamma())
model_gamma.fit(data_gamma, method="REML")
print(model_gamma.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Gamma(link='log')
Observations: 250
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 1.0938 0.0199 55.096 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 6.72 7 178.726 < 1e-16
Total EDF: 7.72
Deviance: 23.8734
Null dev: 41.5510
Dev. expl: 42.5%
GCV score: 0.101680
Scale est: 0.098539
AIC: 665.02
BIC: 692.22
With the log link the model is \log(\mu) = \beta_0 + f(x), so coefficients are interpretable as multiplicative effects on the response.
## Visualizing the Gamma fit
``` python
x_grid = np.linspace(0.5, 5, 300)
preds_gamma = model_gamma.predict({"x": x_grid}, type="response")
points_gamma = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_gamma[i])} for i in range(n)]}
).mark_circle(size=15, opacity=0.3, color="goldenrod").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="y"),
)
line_gamma = alt.Chart(
{"values": [{"x": float(x_grid[i]), "fit": float(preds_gamma.values[i])}
for i in range(len(x_grid))]}
).mark_line(color="firebrick", strokeWidth=2).encode(x="x:Q", y="fit:Q")
(points_gamma + line_gamma).properties(
width="container", height=300,
title="Gamma GAM: fitted mean"
)
```
# Negative Binomial: overdispersed counts
When count data exhibit more variability than the Poisson model allows (i.e., the variance exceeds the mean), the Negative Binomial family provides extra flexibility via an overdispersion parameter \theta:
- **Link**: g(\mu) = \log(\mu).
- **Variance function**: V(\mu) = \mu + \mu^2 / \theta. As \theta \to \infty, this reduces to the Poisson.
- **Deviance**: D = 2\sum_i \bigl\[\theta\\\log\\\bigl(\theta / (\theta + \hat\mu_i)\bigr) + y_i \log\\\bigl(y_i(\theta + \hat\mu_i) / (\hat\mu_i(\theta + y_i))\bigr)\bigr\].
The overdispersion parameter \theta is estimated alongside the smoothing parameters during REML optimization.
**Supported links**: `log` (default), `identity`, `sqrt`.
## Example: overdispersed species counts
``` python
rng = np.random.default_rng(4)
n = 300
x = np.linspace(0, 5, n)
mu = np.exp(1.5 + 0.8 * np.sin(1.5 * x))
theta = 3.0 # overdispersion parameter
# Generate NB data: Poisson-Gamma mixture
lam = rng.gamma(theta, scale=mu / theta, size=n)
y_nb = rng.poisson(lam).astype(float)
data_nb = {"x": x, "y": y_nb}
model_nb = wk.GAM("y ~ s(x, k=15)", family=wk.NegativeBinomial())
model_nb.fit(data_nb, method="REML")
print(model_nb.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x, k=15)
Family: NegativeBinomial(theta=2.919, link='log')
Observations: 300
Coefficients: 15
Parametric coefficients:
Term Estimate Std.Err z value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 1.5588 0.0443 35.215 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x, k=15) 6.46 7 154.264 < 1e-16
Total EDF: 7.46
Deviance: 335.6509
Null dev: 501.0747
Dev. expl: 33.0%
GCV score: 1.176627
Scale est: 1.000000
AIC: 1546.83
BIC: 1574.46
``` python
# Compare observed variance to mean (overdispersion check)
print(f"Observed mean: {y_nb.mean():.2f}")
print(f"Observed variance: {y_nb.var():.2f}")
print(f"Variance / mean: {y_nb.var() / y_nb.mean():.2f}")
```
Observed mean: 5.52
Observed variance: 25.24
Variance / mean: 4.57
A ratio of variance to mean substantially greater than 1 is the hallmark of overdispersion. The Poisson assumption (\text{Var} = \mu) would underestimate standard errors here, leading to overconfident inference. The Negative Binomial family handles this correctly.
## Visualizing the Negative Binomial fit
``` python
x_grid = np.linspace(0, 5, 300)
preds_nb = model_nb.predict({"x": x_grid}, type="response")
points_nb = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_nb[i])} for i in range(n)]}
).mark_circle(size=15, opacity=0.3, color="teal").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="Count"),
)
line_nb = alt.Chart(
{"values": [{"x": float(x_grid[i]), "fit": float(preds_nb.values[i])}
for i in range(len(x_grid))]}
).mark_line(color="darkorange", strokeWidth=2).encode(x="x:Q", y="fit:Q")
(points_nb + line_nb).properties(
width="container", height=300,
title="Negative Binomial GAM: fitted rate"
)
```
# Beta: proportions on (0, 1)
When the response is a continuous proportion (e.g., vegetation cover, exam scores as fractions), the Beta family is the appropriate choice. Unlike the Binomial family (which models counts out of a known total), the Beta family handles continuous values strictly between 0 and 1.
- **Link**: g(\mu) = \log\\\bigl(\mu / (1 - \mu)\bigr) (logit).
- **Variance function**: V(\mu) = \mu(1 - \mu) / (1 + \phi), where \phi is a precision parameter.
- **Deviance**: based on the Beta log-likelihood with parameters \alpha = \mu\phi and \beta = (1-\mu)\phi.
**Supported links**: `logit` (default), `probit`, `cloglog`, `cauchit`.
## Example: proportion of canopy cover
``` python
rng = np.random.default_rng(5)
n = 200
x = np.linspace(0, 4, n)
mu = 1 / (1 + np.exp(-(0.5 + 0.8 * np.sin(1.5 * x)))) # true proportion
phi = 20.0 # precision
alpha = mu * phi
beta_param = (1 - mu) * phi
y_beta = rng.beta(alpha, beta_param)
data_beta = {"x": x, "y": y_beta}
model_beta = wk.GAM("y ~ s(x)", family=wk.Beta())
model_beta.fit(data_beta, method="REML")
print(model_beta.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Beta(link='logit')
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.4953 0.0315 15.724 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 6.14 7 276.307 < 1e-16
Total EDF: 7.14
Deviance: 8.3937
Null dev: 21.1767
Dev. expl: 60.4%
GCV score: 0.045136
Scale est: 0.043523
AIC: -362.73
BIC: -339.16
## Visualizing the Beta fit
``` python
x_grid = np.linspace(0, 4, 300)
preds_beta = model_beta.predict({"x": x_grid}, type="response")
points_beta = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_beta[i])} for i in range(n)]}
).mark_circle(size=15, opacity=0.3, color="seagreen").encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="Proportion"),
)
line_beta = alt.Chart(
{"values": [{"x": float(x_grid[i]), "fit": float(preds_beta.values[i])}
for i in range(len(x_grid))]}
).mark_line(color="darkviolet", strokeWidth=2).encode(x="x:Q", y="fit:Q")
(points_beta + line_beta).properties(
width="container", height=300,
title="Beta GAM: fitted mean proportion"
)
```
> **Warning: Boundary values**
>
> The Beta distribution is defined on the open interval (0, 1). If your data contain exact 0s or 1s, you must transform them slightly (e.g., (y(n-1) + 0.5)/n) or consider a zero/one-inflated model.
# Tweedie: mixed continuous-discrete responses
The Tweedie family generalizes several distributions through a **power parameter** p:
- p = 0: Gaussian
- p = 1: Poisson
- 1 \< p \< 2: compound Poisson-Gamma (a point mass at zero plus a continuous distribution on the positives)
- p = 2: Gamma
- p = 3: inverse Gaussian
The most common use is 1 \< p \< 2, which naturally handles data with exact zeros and a continuous positive tail (e.g., insurance claims, daily rainfall).
- **Link**: g(\mu) = \log(\mu).
- **Variance function**: V(\mu) = \mu^p.
**Supported links**: `log` (default), `identity`.
## Example: insurance claims
``` python
rng = np.random.default_rng(6)
n = 300
x = np.linspace(0, 5, n)
mu = np.exp(0.5 + 0.6 * np.sin(x))
# Simulate Tweedie-like data (compound Poisson-Gamma)
p_tw = 1.5
phi_tw = 1.0
poisson_rate = mu**(2 - p_tw) / ((2 - p_tw) * phi_tw)
gamma_shape = (2 - p_tw) / (p_tw - 1)
gamma_scale = phi_tw * (p_tw - 1) * mu**(p_tw - 1)
n_events = rng.poisson(poisson_rate)
y_tw = np.array([
rng.gamma(gamma_shape, scale=gamma_scale[i], size=n_events[i]).sum()
if n_events[i] > 0 else 0.0
for i in range(n)
])
data_tw = {"x": x, "y": y_tw}
model_tw = wk.GAM("y ~ s(x)", family=wk.Tweedie())
model_tw.fit(data_tw, method="REML")
print(model_tw.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: Tweedie(p=1.5, link='log')
Observations: 300
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 0.6098 0.0502 12.149 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 4.30 5 75.813 6.293e-15
Total EDF: 5.30
Deviance: 295.5233
Null dev: 376.1626
Dev. expl: 21.4%
GCV score: 1.020843
Scale est: 1.002801
AIC: 1007.88
BIC: 1027.52
``` python
# Proportion of exact zeros
print(f"Proportion of zeros: {(y_tw == 0).mean():.2%}")
```
Proportion of zeros: 3.67%
## Visualizing the Tweedie fit
``` python
x_grid = np.linspace(0, 5, 300)
preds_tw = model_tw.predict({"x": x_grid}, type="response")
points_tw = alt.Chart(
{"values": [{"x": float(x[i]), "y": float(y_tw[i]),
"zero": "zero" if y_tw[i] == 0.0 else "positive"}
for i in range(n)]}
).mark_circle(size=15, opacity=0.4).encode(
x=alt.X("x:Q", title="x"),
y=alt.Y("y:Q", title="y"),
color=alt.Color("zero:N", scale=alt.Scale(
domain=["zero", "positive"], range=["crimson", "steelblue"]
), title="Value"),
)
line_tw = alt.Chart(
{"values": [{"x": float(x_grid[i]), "fit": float(preds_tw.values[i])}
for i in range(len(x_grid))]}
).mark_line(color="firebrick", strokeWidth=2).encode(x="x:Q", y="fit:Q")
(points_tw + line_tw).properties(
width="container", height=300,
title="Tweedie GAM: fitted mean with zero mass"
)
```
> **Tip: Estimating the power parameter**
>
> If you are unsure what value of p to use, the [tw()](../reference/tw.md#whittaker.tw) shorthand lets Whittaker estimate the power parameter from the data:
>
> ``` python
> model = wk.GAM("y ~ s(x)", family=wk.tw())
> ```
>
> This adds p to the set of parameters estimated during REML optimization.
# Inverse Gaussian: positive heavy-tailed responses
The Inverse Gaussian family is appropriate for positive continuous data with a heavier right tail than the Gamma. It arises naturally as the distribution of first passage times in Brownian motion.
- **Link**: g(\mu) = 1/\mu^2 (canonical).
- **Variance function**: V(\mu) = \mu^3.
- **Deviance**: D = \sum_i (y_i - \hat\mu_i)^2 / (\hat\mu_i^2 \\ y_i).
**Supported links**: `1/mu^2` (default), `log`, `inverse`, `identity`.
## Example
``` python
from scipy.stats import invgauss as invgauss_dist
rng = np.random.default_rng(7)
n = 200
x = np.linspace(0.5, 4, n)
mu = np.exp(0.5 + 0.3 * x)
y_ig = np.array([invgauss_dist.rvs(mu=m, scale=1.0, random_state=rng) for m in mu])
data_ig = {"x": x, "y": y_ig}
model_ig = wk.GAM("y ~ s(x)", family=wk.InverseGaussian())
model_ig.fit(data_ig, method="REML")
print(model_ig.summary())
```
GAM fit summary
============================================================
Formula: y ~ s(x)
Family: InverseGaussian(link='log')
Observations: 200
Coefficients: 10
Parametric coefficients:
Term Estimate Std.Err t value p-value
------------------------ ---------- ---------- ---------- ----------
(Intercept) 1.2682 0.1334 9.510 < 1e-16
Approximate significance of smooth terms:
Term EDF Ref.df Chi.sq p-value
------------------------ ------ ------ ---------- ----------
s(x) 1.00 2 11.837 0.002689
Total EDF: 2.00
Deviance: 181.1550
Null dev: 193.8270
Dev. expl: 6.5%
GCV score: 0.924166
Scale est: 0.914924
AIC: 807.79
BIC: 814.39
The Inverse Gaussian is less commonly used than the Gamma but should be considered when the coefficient of variation increases more steeply with the mean (since V(\mu) = \mu^3 vs. \mu^2 for Gamma).
# Cox proportional hazards: survival data
The [CoxPH()](../reference/CoxPH.md#whittaker.CoxPH) family supports semiparametric survival analysis via the Cox proportional hazards model. Unlike the other families above, this is not a member of the exponential dispersion family in the traditional sense. Instead, it uses a partial likelihood formulation.
``` python
model = wk.GAM("survival_object ~ s(age) + s(biomarker)", family=wk.CoxPH())
model.fit(data, method="REML")
```
See the full CoxPH documentation in the API reference for details on constructing the survival response object and interpreting hazard ratios.
# Choosing the right family
Selecting a family involves matching the distribution to three characteristics of your response variable:
1. **Support**: what values can y take?
- Unbounded real numbers → Gaussian
- Non-negative integers → Poisson or Negative Binomial
- Binary (0/1) → Binomial
- Strictly positive reals → Gamma, Inverse Gaussian
- Proportions in (0, 1) → Beta
- Non-negative with point mass at zero → Tweedie (1 \< p \< 2)
2. **Mean-variance relationship**: how does the spread change with the level?
- Constant variance → Gaussian
- Variance \propto mean → Poisson
- Variance \propto mean^2 → Gamma
- Variance \> mean (overdispersion) → Negative Binomial
3. **Domain knowledge**: what does the scientific context suggest?
- Rates and counts → Poisson or Negative Binomial
- Times, costs, concentrations → Gamma or Inverse Gaussian
- Survival analysis → CoxPH
> **Tip: Diagnostic checks**
>
> After fitting, use [model diagnostics](diagnostics.md) to verify your family choice. The key checks are:
>
> - **Residual plots**: residuals vs. fitted values should show no systematic pattern.
> - **QQ plot**: quantile-quantile plot of deviance residuals should be approximately linear.
> - **Overdispersion**: for Poisson models, check that the estimated scale parameter is close to 1. If it is much larger, switch to Negative Binomial.
# Where to go next
- **[Smooth terms](smooths.md)**: basis types, tensor products, and choosing `k`.
- **[Model fitting](fitting.md)**: P-IRLS, REML vs. GCV, and convergence settings.
- **[Prediction and inference](prediction.md)**: standard errors, confidence intervals, and term-level predictions.
- **[Diagnostics](diagnostics.md)**: residual plots, QQ plots, and model checking.