Correlation structures

Residual correlation structures for linear mixed models with repeated measures or longitudinal data. Pass a correlation structure to interlace.fit() via the correlation parameter.

The key operation is whitening: transforming (y, X, Z) so the residual covariance becomes identity, allowing the standard profiled-REML machinery to operate on the transformed data. The correlation parameter (rho) is estimated jointly with the variance parameters.

AR(1) — first-order autoregressive

Within each group, residuals separated by time gap dt have correlation rho^dt (continuous-time AR(1) / Ornstein-Uhlenbeck process).

class AR1(time)[source]

First-order autoregressive residual correlation.

Within each group, residuals separated by time gap dt have correlation rho^dt (continuous-time AR(1) / Ornstein-Uhlenbeck process).

Parameters:

time (str) – Name of the column in the data that contains the time/order variable.

Parameters:

time (str)

property n_corr_params: int

Number of correlation parameters to estimate.

setup(groups, times)[source]

Pre-compute group slices, sort order, and time gaps.

Data is sorted by (group, time) internally. The sort and unsort indices are stored so that whitening can be applied efficiently.

Return type:

None

Parameters:
  • groups (ndarray)

  • times (ndarray)

whiten_data(y, X, Z, rho_params)[source]

Apply AR(1) whitening to (y, X, Z).

Parameters:
  • y (shape (n,))

  • X (shape (n, p))

  • Z (sparse (n, q))

  • rho_params (shape (1,)) – The single AR(1) parameter (already mapped from unconstrained space).

Return type:

tuple[ndarray, ndarray, csc_matrix]

Returns:

(y_w, X_w, Z_w) (whitened data in original observation order.)

Parameters:
  • y (ndarray)

  • X (ndarray)

  • Z (csc_matrix)

  • rho_params (ndarray)

log_det_R(rho_params, **kwargs)[source]

Compute log|R(rho)| = sum over groups of sum_{j} log(1 - rho^{2*dt_j}).

When rho=0, log|R| = 0 (R is the identity).

Return type:

float

Parameters:
  • rho_params (ndarray)

  • kwargs (object)

log_det_R_from_arrays(rho, groups, times, sort_idx)[source]

Compute log|R| from raw arrays (used in tests).

Return type:

float

Parameters:
  • rho (float)

  • groups (ndarray)

  • times (ndarray)

  • sort_idx (ndarray)

Example

import interlace
from interlace import AR1

result = interlace.fit(
    formula="y ~ time + treatment",
    data=df,
    groups="subject",
    correlation=AR1(time="time"),
)

# Estimated correlation parameter
print(result.correlation_params)
# {'rho': 0.72}

When to use

AR(1) is appropriate when observations within a group are ordered in time and the correlation decays with increasing temporal distance. Common in longitudinal studies with equally or unequally spaced measurements.


CompoundSymmetry — exchangeable correlation

All pairs of residuals within a group share the same correlation rho. The within-group correlation matrix is R_g = (1-rho)I + rho J, where J is the all-ones matrix.

class CompoundSymmetry(time)[source]

Exchangeable residual correlation within groups.

All pairs of residuals within a group share the same correlation rho. The within-group correlation matrix is R_g = (1-rho)*I + rho*J where J is the all-ones matrix.

Parameters:

time (str) – Name of the grouping column in the data (used by the CorStruct interface; the actual time values are irrelevant for CS).

Parameters:

time (str)

property n_corr_params: int

Number of correlation parameters to estimate.

setup(groups, times)[source]

Pre-compute group slices and sizes.

Data is sorted by group internally so whitening can operate on contiguous slices.

Return type:

None

Parameters:
  • groups (ndarray)

  • times (ndarray)

unconstrained_bounds()[source]

Bounds for the unconstrained rho parameter.

CS requires -1/(n_max - 1) < rho < 1 for positive definiteness. We map the lower limit (with a small margin) to unconstrained space via atanh.

Return type:

list[tuple[float | None, float | None]]

whiten_data(y, X, Z, rho_params)[source]

Apply compound symmetry whitening to (y, X, Z).

Return type:

tuple[ndarray, ndarray, csc_matrix]

Parameters:
  • y (ndarray)

  • X (ndarray)

  • Z (csc_matrix)

  • rho_params (ndarray)

For a group of size n_g with correlation rho:

R^{-1/2} v = v / a - c * sum(v) * 1

where a = sqrt(1 - rho), c = (1/a - 1/sqrt(1 + (n_g-1)*rho)) / n_g.

log_det_R(rho_params, **kwargs)[source]

Compute log|R(rho)| for compound symmetry.

Return type:

float

Parameters:
  • rho_params (ndarray)

  • kwargs (object)

For a group of size n_g:

log|R_g| = (n_g - 1) * log(1 - rho) + log(1 + (n_g - 1) * rho)

Total is the sum over all groups.

Example

from interlace import CompoundSymmetry

result = interlace.fit(
    formula="y ~ time + treatment",
    data=df,
    groups="subject",
    correlation=CompoundSymmetry(time="time"),
)

When to use

Compound symmetry (exchangeable correlation) is appropriate when the within-group correlation does not depend on temporal distance — e.g. when repeated measures are interchangeable or when a random intercept alone does not adequately capture the residual dependence.


Constructor parameters

Both AR1 and CompoundSymmetry accept a single parameter:

Parameter

Type

Description

time

str

Name of the column in the data containing the time/order variable. For AR(1), used to compute time gaps. For CS, used only for the CorStruct interface.

Estimated parameters

After fitting, result.correlation_params is a dict containing:

Key

Description

rho

Estimated within-group correlation parameter

Comparison with R

interlace

R (nlme)

fit(..., correlation=AR1(time="time"))

lme(..., correlation=corAR1(form=~time|subject))

fit(..., correlation=CompoundSymmetry(time="time"))

lme(..., correlation=corCompSymm(form=~time|subject))

result.correlation_params["rho"]

coef(fit$modelStruct$corStruct, unconstrained=FALSE)

See also

  • fit — the correlation parameter on fit()

  • CrossedLMEResultCrossedLMEResult attributes including correlation_params