MultiResponseGAM(
responses,
formula,
*,
response_formulas=None,
family=None,
correlation="independent",
)
Fits multiple response variables jointly against a common set of covariates, optionally sharing smooth terms (same basis and formula structure, though each response still gets its own coefficients) and estimating the residual correlation between responses. Internally, each response is fit as its own GAM using the shared formula plus any response-specific additional terms; what makes this a genuinely multivariate model rather than just several independent fits is the optional joint residual covariance structure, which is useful for understanding how responses co-vary after accounting for the shared covariates, and for joint (GLS-style) prediction via joint_predict().
Use MultiResponseGAM when you have several related outcomes measured on the same units (e.g. multiple biomarkers, or several pollutant concentrations) that likely share similar covariate relationships and whose residuals may be correlated.
Parameters
responses: list[str]-
List of response variable names (at least two).
formula: str-
Shared formula applied to all responses (e.g. "s(x1) + s(x2)"). The response side (before ~, if present) is ignored; use responses to specify the response variables.
response_formulas: dict[str, str] | None = None-
Dict mapping response name to a response-specific formula string (covariates only, e.g., {"y1": "s(x3)"}) added on top of the shared formula for that response only.
family: Family | None = None-
Response distribution family (applied to all responses). Defaults to Gaussian().
correlation: str = "independent"-
Residual correlation structure:
"independent" (default), which fits each response’s GAM independently with no cross-response covariance modeling, or "unstructured", which additionally estimates a full k x k residual covariance matrix from the fitted residuals.
Notes
Under correlation="unstructured", after each response’s GAM is fit, the residual matrix R \in \mathbb{R}^{n \times k} (columns are y_j - \hat y_j for each response j) is used to estimate the residual covariance,
\hat\Sigma = \frac{R^\top R}{n - 1},
and the corresponding correlation matrix by rescaling to unit diagonal. This does not feed back into how the individual response GAMs are fit (each is still fit marginally), but it is used by joint_predict() to report a joint covariance alongside the stacked mean predictions, and by residual_correlation() for diagnosing cross-response dependence.
Examples
import numpy as np
from whittaker.multi_response import MultiResponseGAM
rng = np.random.default_rng(0)
n = 400
x = rng.uniform(0, 1, n)
shared = np.sin(2 * np.pi * x)
y1 = shared + rng.normal(scale=0.2, size=n)
y2 = 0.5 * shared + rng.normal(scale=0.2, size=n)
model = MultiResponseGAM(["y1", "y2"], "s(x)", correlation="unstructured")
model.fit({"x": x, "y1": y1, "y2": y2})
print(model.residual_correlation())
ResidualCorrelation:
y1 y2
y1 1.000 -0.019
y2 -0.019 1.000
Attributes
|
Name
|
Description
|
|
correlation
|
Residual correlation structure used when fitting.
|
|
is_fitted
|
Whether fit() has been called successfully.
|
|
n_responses
|
Number of response variables modeled jointly.
|
|
responses
|
List of response variable names, in the order passed to the constructor.
|
correlation
Residual correlation structure used when fitting.
is_fitted
Whether fit() has been called successfully.
n_responses
Number of response variables modeled jointly.
responses
List of response variable names, in the order passed to the constructor.
Methods
deviance()
Return deviance per response.
edf()
Return total EDF per response.
fit()
Fit the multi-response GAM.
fit(
data,
*,
method="REML",
select=False,
)
Fits one GAM per response, using the shared formula plus that response’s entry (if any) in response_formulas. If correlation="unstructured", also estimates the residual covariance and correlation matrices across responses from the fitted residuals.
Parameters
data: InputData-
Column-oriented data containing all response and covariate columns.
method: str = "REML"-
Smoothing parameter selection method applied to every response’s GAM.
select: bool = False-
Enable double-penalty variable selection for every response’s GAM.
Returns
MultiResponseGAM-
Returns
self for method chaining.
joint_predict()
Predict all responses jointly, returning a matrix.
Parameters
new_data: InputData-
Column-oriented covariate data.
Returns
tuple[NDArray, NDArray | None]-
(predictions, covariance) where predictions is (n, k) and covariance is the (k, k) residual covariance matrix (None if correlation="independent").
predict()
Predict all responses on new data.
predict(
new_data,
*,
se=False,
)
Predicts each response’s fitted GAM independently on the same covariate data, packaging the results into a single MultiResponseResult.
Parameters
new_data: InputData-
Column-oriented covariate data.
se: bool = False-
If
True, include standard errors for each response’s linear predictor.
residual_correlation()
Return the estimated residual correlation structure.
Only available when correlation="unstructured".
response_model()
Return the fitted GAM for a specific response.
Parameters
response: str-
Response variable name.
summary()
Build a text summary of the fitted multi-response GAM.
Reports the response names, shared formula, family, correlation structure, and per-response fit statistics (EDF, deviance, scale). When correlation="unstructured", also lists the estimated pairwise residual correlations between responses.
Returns
str-
Multi-line summary text.