# Multi-response models

When multiple outcomes are measured on the same observations (for example, blood pressure and cholesterol, or multiple pollutant concentrations at the same monitoring stations), fitting them jointly can capture shared structure and estimate residual correlations. Whittaker's [MultiResponseGAM](../reference/MultiResponseGAM.md#whittaker.MultiResponseGAM) fits a separate GAM per response but provides a unified interface for prediction, diagnostics, and correlation estimation.


# Basic multi-response model


``` python
import numpy as np
import whittaker as wk

# Generate two correlated responses driven by the same smooth
rng = np.random.default_rng(23)
n = 300
x = np.linspace(0, 2 * np.pi, n)

e1 = rng.normal(0, 0.3, n)
e2 = 0.7 * e1 + rng.normal(0, 0.2, n)  # correlated errors

y1 = np.sin(x) + e1
y2 = 0.5 * np.sin(x) + 1.5 + e2

data = {"x": x, "y1": y1, "y2": y2}
```


The two responses share the same smooth predictor `x` but have different signal strengths and correlated residuals. We fit them jointly:


``` python
# Fit a multi-response GAM with a shared smooth formula
model = wk.MultiResponseGAM(["y1", "y2"], "s(x)")
model.fit(data, method="REML")

print(f"Responses: {model.responses}")
print(f"N responses: {model.n_responses}")
print(f"Is fitted: {model.is_fitted}")
```


    Responses: ['y1', 'y2']
    N responses: 2
    Is fitted: True


The first argument lists the response column names. The second is the shared formula (covariates only, no response side). Each response gets its own GAM with the same smooth structure.


# Predictions


``` python
# Predict both responses on new data
x_new = np.linspace(0, 2 * np.pi, 100)
result = model.predict({"x": x_new})

# Access predictions by response name
print(f"y1 predictions shape: {result['y1'].values.shape}")
print(f"y2 predictions shape: {result['y2'].values.shape}")
```


    y1 predictions shape: (100,)
    y2 predictions shape: (100,)


``` python
import altair as alt

# Plot both response fits
plot_data = []
for i in range(len(x_new)):
    plot_data.append({"x": float(x_new[i]), "y": float(result["y1"].values[i]), "response": "y1"})
    plot_data.append({"x": float(x_new[i]), "y": float(result["y2"].values[i]), "response": "y2"})

alt.Chart({"values": plot_data}).mark_line(strokeWidth=2).encode(
    x=alt.X("x:Q"),
    y=alt.Y("y:Q"),
    color=alt.Color("response:N"),
).properties(width="container", height=300, title="Multi-response GAM predictions")
```


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

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


# Predictions with standard errors


``` python
# Standard errors for each response
result_se = model.predict({"x": x_new}, se=True)

print(f"y1 SE mean: {result_se['y1'].se.mean():.4f}")
print(f"y2 SE mean: {result_se['y2'].se.mean():.4f}")
```


    y1 SE mean: 0.0520
    y2 SE mean: 0.0463


# Residual correlation

When responses are modeled jointly, the residuals often carry shared information that the smooth terms do not capture. Estimating this residual correlation is useful for understanding the unexplained association between outcomes.


``` python
# Fit with unstructured residual correlation
model_corr = wk.MultiResponseGAM(
    ["y1", "y2"], "s(x)",
    correlation="unstructured",
)
model_corr.fit(data, method="REML")

# Estimate residual correlation
rc = model_corr.residual_correlation()
print(rc)
```


    ResidualCorrelation:
                   y1        y2
          y1    1.000     0.745
          y2    0.745     1.000


The residual correlation matrix shows how much the residuals of each response pair co-vary after accounting for the smooth terms. A strong positive correlation means the unexplained variation in one response tends to go in the same direction as the other.


``` python
# Access the raw matrices
print(f"Covariance matrix:\n{rc.covariance.round(4)}")
print(f"\nCorrelation matrix:\n{rc.correlation.round(4)}")
```


    Covariance matrix:
    [[0.0946 0.0672]
     [0.0672 0.086 ]]

    Correlation matrix:
    [[1.     0.7448]
     [0.7448 1.    ]]


``` python
# Heatmap of the residual correlation matrix
responses = model_corr.responses
records = [
    {"row": r, "col": c, "corr": float(rc.correlation[i, j])}
    for i, r in enumerate(responses)
    for j, c in enumerate(responses)
]

base = alt.Chart({"values": records})

heatmap = base.mark_rect().encode(
    x=alt.X("col:N", title=None),
    y=alt.Y("row:N", title=None),
    color=alt.Color(
        "corr:Q",
        scale=alt.Scale(scheme="redblue", domain=[-1, 1]),
        title="Correlation",
    ),
)

labels = base.mark_text(fontSize=12).encode(
    x="col:N", y="row:N",
    text=alt.Text("corr:Q", format=".3f"),
    color=alt.condition(
        alt.datum.corr > 0.5, alt.value("white"), alt.value("black")
    ),
)

(heatmap + labels).properties(width="container", height=250, title="Residual correlation matrix")
```


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

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


# Joint prediction

For applications that need all responses in a single matrix (e.g., multivariate downstream analysis), use [joint_predict()](../reference/MultiResponseGAM.md#whittaker.MultiResponseGAM.joint_predict):


``` python
# Joint prediction returns (n x k) matrix and optional covariance
preds_matrix, cov_matrix = model_corr.joint_predict({"x": x_new})

print(f"Joint predictions shape: {preds_matrix.shape}")
print(f"Covariance matrix shape: {cov_matrix.shape}")
```


    Joint predictions shape: (100, 2)
    Covariance matrix shape: (2, 2)


# Per-response model access

You can extract the fitted GAM for any individual response for further inspection:


``` python
# Get the individual GAM for y1
gam_y1 = model.response_model("y1")
print(f"y1 GAM fitted: {gam_y1.is_fitted}")
print(f"y1 EDF: {gam_y1.edf_total:.1f}")
```


    y1 GAM fitted: True
    y1 EDF: 8.5


# EDF and deviance per response


``` python
# Effective degrees of freedom per response
edfs = model.edf()
print(f"EDF: {edfs}")

# Deviance per response
devs = model.deviance()
print(f"Deviance: {devs}")
```


    EDF: {'y1': 8.471517794702452, 'y2': 7.42980770486148}
    Deviance: {'y1': 28.282806290464258, 'y2': 25.70879276268304}


# Response-specific formulas

Sometimes different responses need different covariates. Use `response_formulas` to add response-specific terms on top of the shared formula:


``` python
# Add an extra covariate z that only affects y1
rng = np.random.default_rng(23)
z = rng.uniform(0, 1, n)
data_extra = {"x": x, "z": z, "y1": y1 + 2 * z, "y2": y2}

model_specific = wk.MultiResponseGAM(
    ["y1", "y2"],
    "s(x)",
    response_formulas={"y1": "s(z)"},  # y1 gets an extra smooth for z
)
model_specific.fit(data_extra, method="REML")

# y1 should have higher EDF because it has an extra smooth
print(f"y1 EDF: {model_specific.edf()['y1']:.1f}")
print(f"y2 EDF: {model_specific.edf()['y2']:.1f}")
```


    y1 EDF: 9.5
    y2 EDF: 7.4


# Three or more responses

[MultiResponseGAM](../reference/MultiResponseGAM.md#whittaker.MultiResponseGAM) works with any number of responses (minimum 2):


``` python
# Three responses
y3 = np.cos(x) + rng.normal(0, 0.3, n)
data_three = {"x": x, "y1": y1, "y2": y2, "y3": y3}

model_three = wk.MultiResponseGAM(
    ["y1", "y2", "y3"], "s(x)",
    correlation="unstructured",
)
model_three.fit(data_three, method="REML")

# 3x3 correlation matrix
rc3 = model_three.residual_correlation()
print(rc3)
```


    ResidualCorrelation:
                   y1        y2        y3
          y1    1.000     0.745    -0.024
          y2    0.745     1.000     0.001
          y3   -0.024     0.001     1.000


# Summary

The `summary()` method reports per-response fit statistics, including EDF, deviance, and smoothing parameters for each response.


``` python
print(model_corr.summary())
```


    MultiResponseGAM summary
    ============================================================
    Responses:   y1, y2
    Shared:      s(x)
    Family:      Gaussian
    Correlation: unstructured

    Per-response fits:
      y1: edf=8.5, dev=28.3, scale=0.0970
      y2: edf=7.4, dev=25.7, scale=0.0879

    Residual correlations:
      corr(y1, y2) = 0.745


> **Tip: When to use MultiResponseGAM**
>
> Use [MultiResponseGAM](../reference/MultiResponseGAM.md#whittaker.MultiResponseGAM) when:
>
> - You want to fit multiple outcomes with shared covariates in a single call
> - You need to estimate residual correlations between outcomes
> - You want a unified prediction interface for multiple responses
> - Different responses may need different additional terms via `response_formulas`
>
> If the responses are truly independent and you do not need correlation estimates, fitting separate [GAM](../reference/GAM.md#whittaker.GAM) objects is equivalent and may be simpler.

You can now fit multi-response GAMs with shared or response-specific formulas, estimate residual correlations, and access per-response models for further inspection.


# Where to go next

- **[Prediction and inference](prediction.md)**: confidence intervals and term-level predictions for individual response models.
- **[Model diagnostics](diagnostics.md)**: check each response model's residuals and basis adequacy.
- **[Functional regression](functional.md)**: another multi-predictor setting where the covariates are entire curves.