Random Effects

Z Model

Mixed-effects design matrix formulation for defining custom random effect structures using an explicit design matrix Z.

← Back to Random Effects

Overview

The Z model implements the classical mixed-effects design matrix formulation, providing a flexible way to define custom random effect structures. Instead of specifying a grouping variable, you provide an explicit design matrix Z that maps random effects to observations.

This is particularly useful when:

  • You need to implement random effects that don't fit standard models
  • You want to reproduce results from classical mixed-effects software
  • You need fine-grained control over the random effect structure

Mathematical Formulation

The linear predictor in a mixed-effects model can be written as:

\[\eta = X\beta + Z\mathbf{z}\]

where:

  • \(X\) is the \(n \times p\) fixed effects design matrix
  • \(\beta\) is the \(p\)-vector of fixed effect coefficients
  • \(Z\) is the \(n \times m\) random effects design matrix
  • \(\mathbf{z}\) is the \(m\)-vector of random effects

The random effects are assumed to follow:

\[\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \tau^{-1} C^{-1})\]

where \(\tau\) is the precision hyperparameter and \(C\) is a known \(m \times m\) precision structure matrix (default: identity matrix).

Augmented Latent Field

Internally, INLA works with an augmented latent field \(\tilde{\mathbf{z}} = (\mathbf{v}, \mathbf{z})^\top\) of length \(n + m\), where:

  • \(\mathbf{v}\) are auxiliary variables (first \(n\) components)
  • \(\mathbf{z}\) are the actual random effects (last \(m\) components)

The precision structure for this augmented field is:

\[\mathbf{Q} = \kappa \begin{pmatrix} I_n & -Z \\ -Z^\top & Z^\top Z \end{pmatrix} + \tau \begin{pmatrix} 0 & 0 \\ 0 & C \end{pmatrix}\]

where \(\kappa\) is a large precision parameter (default: 108) that ties \(\mathbf{v}\) to \(Z\mathbf{z}\).

When to Use Z Model

The Z model is appropriate when:

Custom Grouping Structures

When observations belong to multiple groups simultaneously (crossed random effects) or have complex nested structures.

Reproducing Classical Results

When matching results from lme4, nlme, or other mixed-effects software that uses explicit Z matrices.

Weighted Random Effects

When random effects should contribute different amounts to different observations (non-binary design matrix entries).

Sparse Structures

When the design matrix is sparse and you want to exploit this for computational efficiency.

Sparse Matrix Support

For large-scale problems, use scipy.sparse matrices for both Z and Cmatrix. This significantly reduces memory usage and improves performance when matrices are sparse:

import scipy.sparse as sp

# Create sparse identity matrix
Z = sp.identity(n, format='csr')

# Or sparse indicator matrix for groups
rows = np.arange(n)
cols = group_ids
data = np.ones(n)
Z = sp.csr_matrix((data, (rows, cols)), shape=(n, n_groups))

Specification in pyINLA

Allowed Parameters

Parameter Required Description
id Yes Column name for observation indices (used for internal bookkeeping)
model Yes Must be 'z'
Z Yes The \(n \times m\) design matrix (numpy array or scipy.sparse matrix)
Cmatrix No The \(m \times m\) precision structure matrix for \(\mathbf{z}\) (numpy array or scipy.sparse matrix). Default: identity matrix.
precision No The \(\kappa\) parameter in the A matrix. Default: 108.
constr No Sum-to-zero constraint on \(\mathbf{z}\). Default: False.
extraconstr No Dict with keys 'A' and 'e' for custom linear constraints \(A\mathbf{z} = e\).
hyper No Prior specification for the precision \(\tau\)

Basic Syntax

# Create design matrix Z
import numpy as np

# Example: 100 observations, 10 groups
n_obs = 100
n_groups = 10
groups = np.repeat(np.arange(n_groups), n_obs // n_groups)

# Build design matrix: Z[i,j] = 1 if observation i in group j
Z = np.zeros((n_obs, n_groups))
for i in range(n_obs):
    Z[i, groups[i]] = 1

# Model specification (intercept included by default)
formula = {
    'response': 'y',
    'fixed': ['x'],
    'random': [
        {
            'id': 'idx',
            'model': 'z',
            'Z': Z
        }
    ]
}

result = pyinla(formula=formula, family='gaussian', data=df)

Relationship to IID Model

The Z model with an identity design matrix \(Z = I_n\) is equivalent to the IID model. This can be useful for verification:

# Z model with identity matrix = IID model
Z_identity = np.eye(n)

model_z = {
    'response': 'y',
    'fixed': ['x'],
    'random': [{'id': 'idx', 'model': 'z', 'Z': Z_identity}]
}

# Equivalent to:
model_iid = {
    'response': 'y',
    'fixed': ['x'],
    'random': [{'id': 'idx', 'model': 'iid'}]
}

Similarly, a grouped Z matrix (indicator matrix where \(Z_{ij} = 1\) if observation \(i\) belongs to group \(j\)) produces results equivalent to an IID model on the group variable.

Custom Precision Structure

The Cmatrix parameter allows you to specify a custom precision structure for the random effects:

# Custom C matrix: correlation between groups
C = np.eye(n_groups)
# Add off-diagonal elements for correlation
for i in range(n_groups - 1):
    C[i, i+1] = 0.3
    C[i+1, i] = 0.3

formula = {
    'response': 'y',
    'fixed': ['x'],
    'random': [{
        'id': 'idx',
        'model': 'z',
        'Z': Z,
        'Cmatrix': C
    }]
}

Note on C Matrix

The C matrix must be positive definite. It represents the precision structure (not covariance), so larger values indicate stronger shrinkage toward zero.

Hyperparameters

The Z model has one hyperparameter:

Parameter Internal Scale Description
\(\tau\) \(\log(\tau)\) Precision of the random effects \(\mathbf{z}\)

Prior Specification

'random': [{
    'id': 'idx',
    'model': 'z',
    'Z': Z,
    'hyper': [{
        'prior': 'pc.prec',      # PC prior for precision
        'param': [1, 0.01],     # P(sigma > 1) = 0.01
        'initial': 4            # Starting value for log(tau)
    }]
}]

Constraints

The Z model supports two types of constraints on the random effects \(\mathbf{z}\):

Sum-to-Zero Constraint

Use constr=True to impose \(\sum_{j=1}^m z_j = 0\):

'random': [{
    'id': 'idx',
    'model': 'z',
    'Z': Z,
    'constr': True  # Sum of z equals zero
}]

Custom Linear Constraints

Use extraconstr to impose custom linear constraints \(A\tilde{\mathbf{z}} = e\) on the augmented latent field:

# Two constraints on augmented field (length n+m)
n, m = Z.shape
A_constr = np.zeros((2, n + m))
A_constr[0, n:n+3] = 1     # Sum of first 3 z elements = 0
A_constr[1, n+3] = 1; A_constr[1, n+4] = -1  # z[4] = z[5]

'random': [{
    'id': 'idx',
    'model': 'z',
    'Z': Z,
    'extraconstr': {
        'A': A_constr,
        'e': [0, 0]  # Both constraints equal zero
    }
}]

Note on Augmented Field

The extraconstr constraint matrix A operates on the augmented latent field of length \(n + m\). The first \(n\) elements are auxiliary variables; the last \(m\) elements are the random effects \(\mathbf{z}\).

Reference

The Z model implementation follows section 3.3.3 "Mixed-effects design matrix (z model)" from:

Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B, 71(2), 319-392.