Model Skewed Positive Data with Gamma

Fit a GAM when the response is a positive continuous variable with right skew.

Use the Gamma family whenever your response is strictly positive and right-skewed (income, insurance claims, reaction times, or wages are classic examples). The log link means each covariate has a multiplicative effect on the expected response. Unlike a log-transformed Gaussian, the Gamma model preserves the original scale and gives well-calibrated prediction intervals. The smooth terms s(age) and s(experience) let the data determine the shape of each relationship.

Load data

Load the built-in wages dataset and inspect its columns.

import whittaker as wk

data = wk.load_dataset("wages", as_frame=True)
data.columns.tolist()

['age', 'experience', 'wage']

The wages dataset has three columns: age, experience, and wage.

data.head()

	age	experience	wage
0	62.323637	9.216817	80.231223
1	42.032395	2.487841	40.024189
2	63.883454	7.251912	60.292852
3	21.799293	3.799293	20.232940
4	46.545724	1.738119	62.383214

Look at the distribution of wage. It is bounded below at zero and has a long right tail, which is exactly the shape the Gamma family is designed for.

Fit

Fit a Gamma GAM with REML smoothness selection. REML is generally preferred for Gamma models because it is less prone to under-smoothing than GCV.

# Fit Gamma GAM with REML selection
model = wk.GAM(
    "wage ~ s(age) + s(experience)",
    family=wk.Gamma(),
).fit(data, method="REML")

model.summary()

GAM fit summary
============================================================
Formula:    wage ~ s(age) + s(experience)
Family:     Gamma(link='log')
Inference:  REML
Observations: 800
Coefficients: 19

Parametric coefficients:
  Term                       Estimate    Std.Err    t value    p-value
  ------------------------ ---------- ---------- ---------- ----------
  (Intercept)                  3.7108     0.0084    440.791    < 1e-16

Approximate significance of smooth terms:
  Term                        EDF Ref.df     Chi.sq    p-value
  ------------------------ ------ ------ ---------- ----------
  s(age)                     4.41      5   2481.673    < 1e-16
  s(experience)              3.44      4    176.897    < 1e-16

Total EDF:  8.85
Scale est:  0.056696
Deviance:   44.8550
Null dev:   205.9584
Dev. expl:  78.2%
GCV score:  0.057331
AIC:        5883.06
BIC:        5924.54

Check the deviance explained and the EDF for each smooth. A high EDF for s(experience) would indicate a strongly non-linear wage-experience profile.

Partial effects

wk.partial_effects(model)

Each panel shows how wage changes with that covariate on the log (link) scale, holding the other covariate at its mean. The shaded band is the 95% confidence interval. A roughly linear partial effect means a log-linear relationship; curvature indicates genuine non-linearity.

Diagnostics

wk.check(model)

For a well-fitting Gamma model, the Q-Q plot should follow the diagonal and the residual-vs-fitted plot should show no obvious funnel pattern.

Predict at new values

Construct a small grid of new covariate values and obtain response-scale predictions with confidence intervals.

import pandas as pd

# Build new-data over age-experience pairs
new_data = pd.DataFrame({
    "age": [25, 35, 45],
    "experience": [2, 10, 20],
})

# Predict wages with confidence intervals
preds = model.predict(new_data, interval="confidence")

# Assemble results with prediction bounds
pd.DataFrame({
    "age": new_data["age"],
    "experience": new_data["experience"],
    "predicted_wage": preds.values,
    "lower": preds.lower,
    "upper": preds.upper,
})

	age	experience	predicted_wage	lower	upper
0	25	2	19.962330	19.201127	20.753709
1	35	10	39.925695	38.243158	41.682256
2	45	20	64.009630	59.996830	68.290819

Predictions are automatically back-transformed from the log scale to the original wage scale, so they are directly interpretable as expected wages in the same units as the response.