Your first distributional model#

Open in Colab Download the notebook

%pip install -q superglm

SuperLSS fits several parameters of a response distribution together. A Gaussian model, for example, can let both the mean and the standard deviation vary with the input columns. Each parameter gets its own predictor.

This walkthrough creates sample data, fits two models and compares their predictions on held-out rows. It needs NumPy, pandas and SuperGLM, with no data download.

Create some data#

Here the response mean depends on age and region. Its standard deviation also increases with age. These are simulated continuous measurements, so a Gaussian response is appropriate.

import numpy as np
import pandas as pd

from superglm import GaussianLS, SuperLSS, cat, s

rng = np.random.default_rng(42)
n = 800
frame = pd.DataFrame(
    {
        "age": rng.uniform(18, 80, n),
        "region": rng.choice(["North", "South"], n),
    }
)
age = frame["age"].to_numpy()
south = (frame["region"] == "South").to_numpy()
mean = 10 + 0.08 * (age - 45) + 0.002 * (age - 45) ** 2 + 1.5 * south
sd = 0.8 + 0.018 * (age - 18)
y = rng.normal(mean, sd)

X_train, X_test = frame.iloc[:600], frame.iloc[600:]
y_train, y_test = y[:600], y[600:]

The split happens before fitting. Both models below learn their spline bases and smoothing parameters from the training rows.

Declare the predictors#

family = GaussianLS()

model = SuperLSS(
    family,
    family.location(s("age", kind="cr", k=8), cat("region")),
    family.scale(s("age", kind="cr", k=6)),
)

The family comes first. Its two helper calls describe the predictors:

  • location(...) models the Gaussian conditional mean. Age has a smooth effect and region has a categorical effect.

  • scale(...) models the Gaussian standard deviation. It has its own smooth age effect and no region effect.

These declarations are configuration. model owns the fit. The age smooths have separate coefficients because they belong to different predictors.

kind="cr" selects a cubic regression spline. k sets its basis size; smoothing determines how much of that flexibility the fit uses. Use a bare string, such as "age", when you want a numeric linear term instead. Categories are always explicit with cat(...), including categories stored as numbers.

Every parameter must have a declaration. An empty call such as family.scale() estimates an intercept-only predictor. It does not fix the parameter to a numeric value, and omitting the call is an error. Each predictor includes an intercept unless you pass intercept=False to its helper.

Use helpers on the same family instance that you pass to SuperLSS. You may put the declarations in any order; names identify their parameters.

Fit the model#

model.fit_reml(X_train, y_train, outer="efs+newton")
print(model.summary())
print(model.diagnose())

  parameter         term       edf         lambda     statistic      rank  \
0  location  (intercept)  1.000000            NaN  24901.056285  1.000000   
1  location          age  4.013779    3666.047317    829.964539  4.857741   
2  location       region  1.000000            NaN    214.113010  1.000000   
3     scale  (intercept)  1.000000            NaN    104.481275  1.000000   
4     scale          age  2.299354  126553.268741     81.534083  2.828469   

         p_value   estimate        se note  
0   0.000000e+00  12.377733  0.078439       
1  6.241621e-177        NaN       NaN       
2   1.740059e-48  -1.513514  0.103434       
3   1.586860e-24   0.297882  0.029142       
4   5.280243e-18        NaN       NaN       
SuperLSS fit diagnosis — family: GaussianLS; status: converged_certified; revision: 1

Rows: 600    Coefficients: 15    Fit time: 122.0 ms

Work
Outer EFS iterations       4
Coefficient fits           5
Total inner iterations     16
Ordinary outer proposals rejected   0
Ordinary outer proposals backtracked 0

Time distribution
Phase                         Time    Share   Calls
layout penalty assembly        79.9 ms   65.5%       6
orchestration and unmeasured   14.4 ms   11.8%       0
efs update backtracking         7.3 ms    6.0%      14
predictor compilation           5.8 ms    4.8%       1
likelihood evaluation           5.5 ms    4.5%      33
coefficient decomposition solve   5.1 ms    4.2%      16
curvature gradient assembly     2.0 ms    1.6%      17
terminal inference and null fit   1.9 ms    1.6%       1
dense predictor matrices        134 µs    0.1%       1
frame normalization              11 µs    0.0%       1

Smoothing parameters
Component                         Initial λ    Final λ  Moves  Lead share  Cap iters  Term EDF  Result
location:age#wiggle                4.728e+04       3666      4      25.0%         0    4.014  finite
scale:age#wiggle                   5.400e+05  1.266e+05      4      75.0%         0    2.299  finite

No solver pathology was detected in the available evidence.

Scope: retained fit telemetry only; use detail='full' for evidence, caveats and limitations.

fit_reml estimates the coefficients and smoothing parameters jointly. outer="efs+newton" adds Newton refinement to the smoothing updates. On this example, the default EFS updates stop after rejecting a proposal; the refinement reaches a stationary fit with the same convergence tolerances. diagnose() reports the stopping evidence. Use fit when you want to hold smoothing parameters fixed. Both fitting methods update the model and return it.

The summary separates terms by predictor. An age effect in location changes the conditional mean. An age effect in scale changes the spread. The Gaussian scale link is log(scale - scale_floor), so its coefficients act on that linked quantity.

Predict the mean, spread and a quantile#

parameters = model.predict_parameters(X_test)
predicted_mean = model.predict(X_test)
upper_quantile = model.predict_quantile(X_test, 0.95)

predictions = parameters.assign(
    observed=y_test,
    predicted_mean=predicted_mean,
    q95=upper_quantile,
)
print(predictions.head())

      location     scale   observed  predicted_mean        q95
600  10.155272  1.417398   8.819948       10.155272  12.486685
601  12.261352  1.765724  12.195464       12.261352  15.165710
602  13.923016  1.826905  11.803275       13.923016  16.928007
603  11.863718  1.733140  13.391061       11.863718  14.714479
604  14.724513  1.835942  13.037870       14.724513  17.744369

For this family, parameters has location and scale columns. They contain the mean and standard deviation on the response scale. predict returns the same mean as the location column. predict_link is available when you need the linear predictors before applying the inverse links.

The 95th percentile describes the upper part of each row’s predictive distribution. It is not a confidence bound on the estimated mean.

Compare against constant spread#

Keep the same mean terms and estimate one standard deviation for all rows:

constant_scale = SuperLSS(
    family,
    family.location(s("age", kind="cr", k=8), cat("region")),
    family.scale(),
).fit_reml(X_train, y_train, outer="efs+newton")

held_out_loss = pd.Series(
    {
        "varying_scale": model.scores(X_test, y_test, which=("log",))["log"].mean(),
        "constant_scale": constant_scale.scores(X_test, y_test, which=("log",))["log"].mean(),
    },
    name="mean_negative_log_likelihood",
)
print(held_out_loss)

varying_scale     1.706871
constant_scale    1.785341
Name: mean_negative_log_likelihood, dtype: float64

A lower mean negative log-likelihood is better on these held-out rows. It scores the predicted distribution, so spread matters as well as mean. For your own data, make the split respect time or group boundaries when those matter.

Reusing family does not share fitted coefficients. Constructing and fitting constant_scale leaves the first model’s fit intact. The models each copy their configuration at construction.

Choose another response family#

Choose the family for the response you have. GammaLS models strictly positive values. TweedieLSS admits both zeros and positive values. NegativeBinomialLS models overdispersed counts. Their predictor names differ because their parameters differ.

For example, a Tweedie declaration uses three helpers:

from superglm import TweedieLSS

tweedie = TweedieLSS()
tweedie_model = SuperLSS(
    tweedie,
    tweedie.mu(s("age", kind="cr", k=8), cat("region")),
    tweedie.phi(s("age", kind="cr", k=6)),
    tweedie.p(),
)

This declares the mean, dispersion and power predictors. Their names in results and offsets are mean, dispersion and power. The empty p() call estimates one power value for all rows. This code only constructs the model; fit it to a response for which the Tweedie law is appropriate.

See family predictor names for all nine families, including what their scale and shape parameters mean. For fitting options and return values, use the SuperLSS API reference. The distributional model guide covers weights, offsets, interactions and discrete fitting; the checking guide covers calibration and predictive diagnostics.