Fit a Multi-Response GAM

Model multiple correlated responses simultaneously with a shared covariate structure.

When multiple outcomes are measured on the same observations and share common predictors, modeling them jointly can improve efficiency and reveal how the responses relate to each other. A multi-response GAM fits a separate smooth model for each outcome while optionally estimating the residual correlation structure between them.

Fitting with correlation="unstructured" goes further: after the smooth effects are removed, the remaining residual covariance is estimated, exposing shared variation that the covariates do not explain.

Simulate data

We generate two measurements along a periodic covariate. Both follow sinusoidal patterns but with different phases.

import numpy as np
import whittaker as wk

# Generate periodic covariate
rng = np.random.default_rng(23)
n = 400
x = np.linspace(0, 2 * np.pi, n)

# Simulate two correlated responses
y1 = np.sin(x) + rng.normal(0, 0.3, n)
y2 = np.cos(x) + 0.5 * np.sin(x) + rng.normal(0, 0.3, n)
data = {"y1": y1, "y2": y2, "x": x}

Fit

The formula argument specifies the covariate structure shared by all responses. Each response gets its own smooth fitted independently.

model = wk.MultiResponseGAM(
    responses=["y1", "y2"],
    formula="s(x, k=10)",
).fit(data)
model.summary()

'MultiResponseGAM summary\n============================================================\nResponses:   y1, y2\nShared:      s(x, k=10)\nFamily:      Gaussian\nCorrelation: independent\n\nPer-response fits:\n  y1: edf=8.6, dev=37.2, scale=0.0950\n  y2: edf=8.9, dev=34.5, scale=0.0882'

Predict

.predict() returns a MultiResponseResult whose predictions are accessed by response name. Each entry holds the fitted values as a plain array.

# Build prediction grid and predict
new_data = {"x": np.linspace(0, 2 * np.pi, 100)}
result = model.predict(new_data)
result["y1"].values[:5]

array([0.00353903, 0.06418029, 0.12469097, 0.18480638, 0.24415738])

result["y2"].values[:5]

array([1.00951036, 1.03605396, 1.06231726, 1.08772635, 1.11146978])

Residual correlation

Refitting with correlation="unstructured" estimates the full residual covariance matrix. .residual_correlation() returns the covariance and the standardized correlation alongside the response names.

# Refit with unstructured correlation
model_corr = wk.MultiResponseGAM(
    responses=["y1", "y2"],
    formula="s(x, k=10)",
    correlation="unstructured",
).fit(data)

# Extract residual correlation matrix
rc = model_corr.residual_correlation()
rc.correlation

array([[ 1.        , -0.03647635],
       [-0.03647635,  1.        ]])

Interpret

The rc.correlation matrix shows the residual correlation between y1 and y2 after the smooth effect of x has been removed. A positive off-diagonal value means the two measurements tend to deviate from their respective smooth curves in the same direction. This indicates that something beyond x drives both outcomes together. This shared driver could represent an unmeasured covariate, individual-level random effects, or genuine coupling between the physiological processes being measured.