RW1 Model

Overview

The RW1 model (Random Walk of order 1) is a fundamental model for capturing smooth temporal or spatial variation. It assumes that:

• First differences are independent - \(x_{i+1} - x_i \sim N(0, 1/\tau)\)
• Adjacent values are correlated - nearby time points or locations share information
• Intrinsically Gaussian - forms a Gaussian Markov Random Field (GMRF)

RW1 is ideal for modeling trends and piecewise linear patterns in ordered data, such as time series, age effects, or spatial transects.

Mathematical Formulation

For \(n\) ordered positions, the RW1 model specifies:

\[x_{i+1} - x_i \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \tau^{-1}), \quad i = 1, \ldots, n-1\]

This can be written equivalently as:

\[x_{i+1} \mid x_i \sim \mathcal{N}(x_i, \tau^{-1})\]

where:

• \(x_i\) is the random effect at position \(i\)
• \(\tau\) is the precision of the increments (higher = smoother)
• \(\sigma = 1/\sqrt{\tau}\) is the standard deviation of the increments

The precision matrix for the full vector \(\mathbf{x} = (x_1, \ldots, x_n)^\top\) is:

\[\mathbf{Q} = \tau \begin{pmatrix} 1 & -1 & & \\ -1 & 2 & -1 & \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 1 \end{pmatrix}\]

This tridiagonal structure makes RW1 computationally efficient, even for large \(n\).

When to Use RW1

The RW1 model is appropriate when you have ordered data where:

Visual Understanding: How RW1 Works

RW1 creates a correlated random walk by linking adjacent values: each step is a Gaussian increment whose variance is controlled by the precision \(\tau\). Higher precision means smaller increments and a tighter trajectory; lower precision allows the walk to wander further between adjacent positions.

Random Walk Realizations

Each curve below is a single realization of an RW1 process. At high precision the trajectory stays tight; at low precision it can wander, with occasional large excursions.

Effect of Precision on Smoothness

The precision \(\tau\) controls how much the curve can wiggle. Higher precision = smoother trends.

First Differences: The Building Block

RW1 is defined by its first differences \(\Delta x_i = x_{i+1} - x_i\). These differences are IID normal - the model's single hyperparameter controls their variance.

RW1 vs IID: Side-by-Side Comparison

Compare the correlated RW1 (each step depends on the previous) with IID (every value drawn independently). The same number of parameters produces very different patterns.

Correlation Structure: Why RW1 has Local Dependence

The key difference between IID and RW1 is their correlation structure. In RW1, nearby positions share information - earlier positions "remember" their influence on later ones.

Linear vs Cyclic: Topology of the Neighbor Graph

Setting 'cyclic': True changes the neighbor graph from a line to a circle: position n becomes a neighbor of position 1. This adds one extra increment, x₁ − x_n, to the prior, so values at the start and end of the index are tied together. Same single hyperparameter τ, same intrinsic structure, just one extra link in the graph.

A sample path of 50 points makes the difference visible: the linear walk drifts (start and end fall at different heights), while the cyclic walk is constrained to come back to its starting value.

When Data Looks Like...

Specification in pyINLA

Allowed Parameters

The RW1 random-effect spec dict accepts the following keys. The source of truth is model_defs/latent/rw1.py's ALLOWED_KEYS set, mirrored 1:1 here.

Per-Parameter Runnable Examples

Click any parameter to expand a minimal runnable snippet. All snippets assume the shared setup block below (\(n = 30\) positions, \(\text{nobs} = 80\) observations on an integer time index 1..30).

# shared setup for every snippet below
import numpy as np
import pandas as pd
from pyinla import pyinla

rng = np.random.default_rng(42)
n    = 30           # number of ordered positions
nobs = 80           # observations
df = pd.DataFrame({
    'y':    rng.standard_normal(nobs),
    'x':    rng.standard_normal(nobs),
    'time': rng.integers(1, n + 1, size=nobs),
})

# Bare rw1: 'id' is the column holding the ordered position index.
# 'n' is auto-detected from max(id); constr=True is on by default.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1'}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Default: prior='loggamma', param=[1, 0.00005], initial=4, fixed=False.
# Below: a more informative loggamma and a different starting value.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'hyper': [{'prior': 'loggamma',
                            'param': [1, 0.01],
                            'initial': 2,
                            'fixed': False}]}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Set 'n' larger than max(id) AND supply 'values' to pin the explicit
# location grid. Latent positions 1..40 are estimated even though only
# positions 1..30 appear in the data.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'n': n + 10,
                 'values': list(range(1, n + 11))}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Default is constr=True. Disabling it leaves the latent field's level
# unidentifiable, so the intercept ('1') MUST be removed from 'fixed'.
formula = {
    'response': 'y',
    'fixed': ['0', 'x'],          # no intercept
    'random': [{'id': 'time', 'model': 'rw1',
                 'constr': False}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# When cyclic=True the latent field wraps: x_n is a neighbor of x_1.
# The sum-to-zero constraint is removed automatically by the engine.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'cyclic': True}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Rescales the precision so the latent field's generalized variance
# is 1 regardless of n. Useful for making the prior on the
# log-precision comparable across rw1 components of different sizes.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'scale.model': True}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Pin the full ordered set of positions explicitly. The list may extend
# beyond the observed range to enable prediction at new positions.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'values': list(range(1, n + 11))}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# The expected number of columns in A depends on how the latent grid is anchored:
#
# (A) Auto-detected grid (no 'values', no 'n'):
#       n_eff = len(unique(time))
#       A has shape (k, n_eff - 1)   # sum-to-zero augmentation reduces dim by 1
#     This is data-dependent: missing positions in 'time' shrink n_eff.
#
# (B) Anchored grid (supply BOTH 'n' AND 'values'):
#       n_eff = n   (length of the values list)
#       A has shape (k, n)           # full field dimension
#     Robust against sparse data; recommended for prediction-style usage.
#
# Below: case (B) pins the first and last positions to zero. With n=30,
# A has shape (2, 30); 'e' has length 2.
A = np.zeros((2, n))
A[0, 0]     = 1.0
A[1, n - 1] = 1.0
e = np.array([0.0, 0.0])
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'n':           n,
                 'values':      list(range(1, n + 1)),
                 'extraconstr': {'A': A, 'e': e}}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Small numerical-stability jitter. Rarely needed for rw1 with a proper prior.
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'diagonal': 1e-4}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

# Per-observation weight vector applied to the rw1 contribution.
# Length must equal the number of observations.
w = rng.uniform(0.5, 1.5, size=nobs)
formula = {
    'response': 'y',
    'fixed': ['1', 'x'],
    'random': [{'id': 'time', 'model': 'rw1',
                 'weights': w}],
}
result = pyinla(formula=formula, family='gaussian', data=df)

Index Handling

The id column should contain integer values representing ordered positions. pyINLA maps these to consecutive indices internally:

Key points:

• Indices define the order - the RW1 structure depends on which positions are adjacent
• Multiple observations can share the same index (e.g., multiple measurements at the same time point)
• Gaps in indices are handled automatically (positions are mapped consecutively)

Missing values (NaN)

pyINLA accepts np.nan / NA in two places without error:

• NaN in the time column: that observation contributes nothing to the rw1 component. The covariates file records a -1 sentinel for those rows. The remaining observations drive the latent field as usual.
• NaN in the weights vector: also accepted. The affected rows contribute no weighted rw1 effect.

Both NaN positions and NaN weights interact cleanly with the other rw1 keys (constr, cyclic, explicit values). See the nan_in_idx_* and nan_in_weights cases in the rw1 parity suite for the verified combinations.

Hyperparameters

The RW1 model has one hyperparameter: the precision \(\tau\) of the first differences. In pyINLA, we work with the log-transformed precision \(\theta = \log(\tau)\).

Default Configuration

Customizing Hyperparameters

Override defaults using the hyper key:

'random': [
    {
        'id': 'time',
        'model': 'rw1',
        'hyper': [{
            'prior': 'loggamma',
            'param': [1, 0.01],      # More informative prior
            'initial': 2,             # Different starting point
            'fixed': False
        }]
    }
]

Available Priors

Log-Gamma Prior (Default)

The default loggamma prior places a gamma distribution on the precision \(\tau\):

\[\tau \sim \text{Gamma}(a, b)\]

With default param=[1, 0.00005], this is a weakly informative prior allowing both rough and smooth curves.

# Default loggamma prior
'hyper': [{'prior': 'loggamma', 'param': [1, 0.00005]}]

# More informative: expect moderate smoothing
'hyper': [{'prior': 'loggamma', 'param': [1, 0.01]}]

PC Prior (Penalised Complexity)

The pc.prec prior penalizes deviation from the base model (constant effect). Parameterized via the marginal standard deviation:

\[P(\sigma > u) = \alpha\]

# PC prior: P(sigma > 1) = 0.01
# "I believe the increment SD is unlikely to exceed 1"
'hyper': [{'prior': 'pc.prec', 'param': [1, 0.01]}]

# PC prior: P(sigma > 0.5) = 0.05
# Expect a smooth, slowly-varying trend
'hyper': [{'prior': 'pc.prec', 'param': [0.5, 0.05]}]

Fixed Precision

To fix the precision at a known value:

import math

# Fix precision at tau = 100 (very smooth)
'hyper': [{
    'initial': math.log(100),  # log(tau) = 4.605
    'fixed': True
}]

Sum-to-Zero Constraint

By default, RW1 includes a sum-to-zero constraint (constr=True) for identifiability with an intercept. You can disable this:

# Disable constraint (remove intercept from model)
model = {
    'response': 'y',
    'fixed': ['0', 'x'],  # No intercept
    'random': [
        {
            'id': 'time',
            'model': 'rw1',
            'constr': False
        }
    ]
}

Interpreting Results

Hyperparameter Summary

After fitting, access the precision estimate via result.summary_hyperpar:

print(result.summary_hyperpar)
#                      mean        sd  0.025quant  0.5quant  0.975quant     mode
# Precision for time  45.23   12.34      24.56     43.21      74.89     41.02

Higher precision = smoother trend (smaller increments). To convert to the posterior of the standard deviation of increments, transform the precision marginal with pyinla.tmarginal. Plugging 1/sqrt(mean(tau)) in directly is biased: by Jensen's inequality the convex map \(\tau \mapsto 1/\sqrt{\tau}\) gives \(1/\sqrt{E[\tau]} \le E[1/\sqrt{\tau}]\), so the naive expression underestimates \(E[\sigma]\).

tau_marg = result.marginals_hyperpar['Precision for time']
sigma_marg = pyinla.tmarginal(lambda tau: 1 / np.sqrt(tau), tau_marg)
sigma_summary = pyinla.zmarginal(sigma_marg)  # mean, sd, quantiles of the increment SD
print(f"Posterior mean of increment SD: {sigma_summary['mean']:.4f}")

For a quick "central value" of \(\sigma\) only, the median is exact under any monotone transform (no marginal transform needed):

tau_median = result.summary_hyperpar['0.5quant'].values[0]
sigma_median = 1 / np.sqrt(tau_median)
print(f"Posterior median of increment SD: {sigma_median:.4f}")

Random Effect Estimates

The posterior estimates for each time point are in result.summary_random:

# Temporal trend estimates
trend = result.summary_random['time']
print(trend)
#    ID      mean        sd  0.025quant  0.5quant  0.975quant      mode
# 0   1 -0.234567  0.089123   -0.409258 -0.234567   -0.059876 -0.234567
# 1   2 -0.123456  0.081234   -0.282345 -0.123456    0.035433 -0.123456
# ...

# Plot the trend
import matplotlib.pyplot as plt
plt.fill_between(trend['ID'], trend['0.025quant'], trend['0.975quant'], alpha=0.3)
plt.plot(trend['ID'], trend['mean'], 'b-')
plt.xlabel('Time')
plt.ylabel('RW1 Effect')