# Fit a Multi-Response GAM

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.


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


``` python
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](../reference/MultiResponseResult.md#whittaker.MultiResponseResult) whose predictions are accessed by response name. Each entry holds the fitted values as a plain array.


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


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


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