# Model Skewed Positive Data with Gamma

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.


``` python
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`.


``` python
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.


``` python
# 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


``` python
wk.partial_effects(model)
```


<style>
  #altair-viz-3391ea64d13c4ae0b42af47a9ac17e1a.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-3391ea64d13c4ae0b42af47a9ac17e1a.vega-embed details,
  #altair-viz-3391ea64d13c4ae0b42af47a9ac17e1a.vega-embed details summary {
    position: relative;
  }
</style>


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


``` python
wk.check(model)
```


<style>
  #altair-viz-a4204cb7189d4ba5ad7923eaacf1681b.vega-embed {
    width: 100%;
    display: flex;
  }

  #altair-viz-a4204cb7189d4ba5ad7923eaacf1681b.vega-embed details,
  #altair-viz-a4204cb7189d4ba5ad7923eaacf1681b.vega-embed details summary {
    position: relative;
  }
</style>


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.


``` python
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.