Overview
The RW2 model (Random Walk of order 2) is a smoothing prior that penalizes deviations from linearity. It assumes that:
- Second differences are independent - \(x_{i-1} - 2x_i + x_{i+1} \sim N(0, 1/\tau)\)
- Penalizes curvature - produces smoother, more gradual changes than RW1
- Intrinsic GMRF - sparse, neighbor-only precision matrix \(\mathbf{Q}\) with two null directions (a constant and a linear trend), so the prior is improper without constraints
RW2 is ideal for modeling smooth curves and gradual trends in ordered data. It acts as a smoothing-spline prior, producing curves that are locally linear with smooth transitions.
Mathematical Formulation
For \(n\) ordered positions, the RW2 model specifies:
\[x_{i-1} - 2x_i + x_{i+1} \stackrel{\text{iid}}{\sim} \mathcal{N}(0, \tau^{-1}), \quad i = 2, \ldots, n-1\]
This penalizes the second derivative (curvature), meaning:
- \(x_i\) is the random effect at position \(i\)
- \(\tau\) is the precision of the second differences (higher = smoother)
- \(\sigma = 1/\sqrt{\tau}\) is the standard deviation of the curvature
The precision matrix for the full vector \(\mathbf{x} = (x_1, \ldots, x_n)^\top\) is:
\[\mathbf{Q} = \tau \begin{pmatrix} 1 & -2 & 1 & & \\ -2 & 5 & -4 & 1 & \\ 1 & -4 & 6 & -4 & 1 \\ & \ddots & \ddots & \ddots & \ddots & \ddots \\ & & 1 & -4 & 5 & -2 \\ & & & 1 & -2 & 1 \end{pmatrix}\]
This pentadiagonal structure makes RW2 computationally efficient while producing smoother curves than RW1.
Identifiability Note
RW2 is an improper prior - the precision matrix is rank-deficient (rank \(n-2\)). The null space includes both a constant and a linear trend, so two directions are unidentified. By default, INLA adds a sum-to-zero constraint (constr=True) which removes the constant null direction. The linear-trend null direction is left unconstrained and is informed by the data through prior shrinkage; if you also want to pin the linear trend (e.g. when a time covariate is in the model), add an extraconstr for it. If you disable constr entirely, also remove the intercept from your model.
When to Use RW2
The RW2 model is appropriate when you expect smooth, gradually changing patterns:
Smooth Temporal Trends
Capture gradual temporal evolution without abrupt jumps. Ideal for long-term trends in climate, epidemiology, or economic data where changes are expected to be smooth.
Seasonal Patterns
Model smooth seasonal effects (with cyclic=True) where the pattern wraps around, such as monthly or weekly cycles.
Dose-Response Curves
Estimate smooth relationships between dose/exposure and response without specifying a functional form. RW2's locally-linear posterior shape can adapt to monotone, plateauing, or unimodal patterns as the data warrant.
Nonparametric Smoothing
Semi-parametric regression with a flexible, smooth function. RW2 produces curves similar to cubic smoothing splines.
Use RW1 when you expect locally constant behavior: a random walk whose level can shift between adjacent positions (penalizes first differences). Use RW2 when you expect locally linear behavior: smoother curves where the slope changes gradually with no abrupt kinks (penalizes second differences). RW2 is the right choice when the underlying process is expected to be continuous and roughly differentiable.
Specification in pyINLA
Allowed Parameters
The RW2 random-effect spec dict accepts the following keys. The source of truth is model_defs/latent/rw2.py's ALLOWED_KEYS set, mirrored 1:1 here.
| Parameter | Required | Description |
|---|---|---|
| Core (model-specific) | ||
model |
Yes | Must be 'rw2'. |
id |
Yes | Column name containing time/position indices (integers). |
hyper |
No | List with one dict: [{prior, param, initial, fixed}] for the log-precision \(\theta\). |
n |
No | Number of positions (auto-detected if omitted). |
constr |
No | Sum-to-zero constraint. Default: True. |
cyclic |
No | Circular/periodic random walk. Default: False. |
scale.model |
No | Scale precision matrix to have generalized variance 1. Default: False. |
values |
No | Custom location grid (list of integers): allows prediction at unobserved locations. |
extraconstr |
No | Extra linear constraints: {'A': matrix, 'e': vector} where \(\mathbf{A}\mathbf{x} = \mathbf{e}\). |
diagonal |
No | Extra constant added to the precision-matrix diagonal for numerical stability. Default: 0.0. |
| Cross-model common keys | ||
weights |
No | Per-observation weighting vector applied to the random-effect contribution. Length must equal the number of observations. |
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),
})id + model required pair (minimal call)
# Bare rw2: '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': 'rw2'}],
}
result = pyinla(formula=formula, family='gaussian', data=df)hyper override the loggamma prior on log-precision
# 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': 'rw2',
'hyper': [{'prior': 'loggamma',
'param': [1, 0.01],
'initial': 2,
'fixed': False}]}],
}
result = pyinla(formula=formula, family='gaussian', data=df)n + values extend the position grid for prediction at unobserved locations
# 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': 'rw2',
'n': n + 10,
'values': list(range(1, n + 11))}],
}
result = pyinla(formula=formula, family='gaussian', data=df)constr disable the sum-to-zero constraint (also drop the intercept)
# Default is constr=True (removes the constant null direction).
# Disabling it leaves the latent field's level unidentifiable,
# so the intercept ('1') MUST be removed from 'fixed'. Note: rw2
# also has a linear-trend null direction that constr does NOT remove;
# add an extraconstr if you also want to pin the linear trend.
formula = {
'response': 'y',
'fixed': ['0', 'x'], # no intercept
'random': [{'id': 'time', 'model': 'rw2',
'constr': False}],
}
result = pyinla(formula=formula, family='gaussian', data=df)cyclic circular random walk (position n wraps to position 1)
# 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': 'rw2',
'cyclic': True}],
}
result = pyinla(formula=formula, family='gaussian', data=df)scale.model scale Q to generalized variance 1
# Rescales the precision so the latent field's generalized variance
# is 1 regardless of n. Strongly recommended when combining rw2
# with PC priors so the prior on log-precision is comparable across
# rw2 components of different sizes.
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [{'id': 'time', 'model': 'rw2',
'scale.model': True}],
}
result = pyinla(formula=formula, family='gaussian', data=df)values explicit location grid
# 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': 'rw2',
'values': list(range(1, n + 11))}],
}
result = pyinla(formula=formula, family='gaussian', data=df)extraconstr additional linear constraints \(\mathbf{A}\mathbf{x} = \mathbf{e}\)
# 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 linear-trend null direction by forcing the
# sum of i * x_i to zero (rw2 leaves this direction unconstrained by default).
# With n=30, A has shape (1, 30); 'e' has length 1.
A = np.arange(1, n + 1, dtype=float).reshape(1, n)
e = np.array([0.0])
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [{'id': 'time', 'model': 'rw2',
'n': n,
'values': list(range(1, n + 1)),
'extraconstr': {'A': A, 'e': e}}],
}
result = pyinla(formula=formula, family='gaussian', data=df)diagonal extra constant on the precision diagonal
# Small numerical-stability jitter. Rarely needed for rw2 with a proper prior.
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [{'id': 'time', 'model': 'rw2',
'diagonal': 1e-4}],
}
result = pyinla(formula=formula, family='gaussian', data=df)weights per-observation weighting
# Per-observation weight vector applied to the rw2 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': 'rw2',
'weights': w}],
}
result = pyinla(formula=formula, family='gaussian', data=df)Hyperparameters
The RW2 model has one hyperparameter: the precision \(\tau\) of the second differences. In pyINLA, we work with the log-transformed precision \(\theta = \log(\tau)\).
Default Configuration
| Parameter | Default Value | Description |
|---|---|---|
prior |
'loggamma' |
Log-gamma prior on precision (gamma prior on \(\tau\)) |
param |
[1, 0.00005] |
Shape and rate parameters for the gamma prior |
initial |
4 |
Starting value for \(\theta = \log(\tau)\); corresponds to \(\tau \approx 54.6\) |
fixed |
False |
Whether to estimate (False) or fix (True) the precision |
Customizing Hyperparameters
Override defaults using the hyper key:
'random': [
{
'id': 'time',
'model': 'rw2',
'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 (linear effect). This is the recommended prior especially when using scale.model=True:
\[P(\sigma > u) = \alpha\]
# PC prior: P(sigma > 1) = 0.01
# "I believe the curvature SD is unlikely to exceed 1"
'hyper': [{'prior': 'pc.prec', 'param': [1, 0.01]}]
# PC prior: P(sigma > 0.5) = 0.05
# Expect a very smooth curve
'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
}]Cyclic (Periodic) RW2
Set cyclic=True for circular data where the last position connects back to the first. This is useful for:
- Monthly seasonal effects (month 12 connects to month 1)
- Weekly patterns (day 7 connects to day 1)
- Circular spatial data (e.g., compass directions)
# Monthly seasonal effect with cyclic RW2
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [
{
'id': 'month',
'model': 'rw2',
'cyclic': True
}
]
}
result = pyinla(formula=formula, family='gaussian', data=df)With cyclic=True, the precision matrix wraps around, connecting position \(n\) to position 1 smoothly.
Scaled Precision Matrix
Set scale.model=True to scale the precision matrix so that the generalized variance equals 1. This is highly recommended when using PC priors, as it makes the prior specification more interpretable.
# Recommended setup: scaled RW2 with PC prior
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [
{
'id': 'time',
'model': 'rw2',
'scale.model': True,
'hyper': [{'prior': 'pc.prec', 'param': [1, 0.01]}]
}
]
}
result = pyinla(formula=formula, family='gaussian', data=df)When using PC priors, always set scale.model=True. This ensures the prior parameters have consistent interpretation regardless of the model size \(n\).
Custom Location Grid
The values parameter allows you to specify a custom grid of locations, which is useful for:
- Prediction - include locations where you want predictions but have no data
- Non-contiguous data - data observed at irregular time points
- Extended grid - model effects beyond the observed data range
# Extend the latent grid to predict beyond observed data.
# Data has t = 1:50, but we want effects for t = 1:100
formula = {
'response': 'y',
'fixed': ['1', 'x'],
'random': [
{
'id': 'time',
'model': 'rw2',
'values': list(range(1, 101)) # [1, 2, ..., 100]
}
]
}
result = pyinla(formula=formula, family='gaussian', data=df)Important: All data values in the id column must be present in the values list. The model will have effects at all specified locations, with data contributing only where observations exist.
Sum-to-Zero Constraint
By default, RW2 includes a sum-to-zero constraint (constr=True) for identifiability with an intercept. You can disable this:
# Disable constraint (remove intercept from formula)
formula = {
'response': 'y',
'fixed': ['0', 'x'], # No intercept
'random': [
{
'id': 'time',
'model': 'rw2',
'constr': False
}
]
}
result = pyinla(formula=formula, family='gaussian', data=df)When to Disable Constraint
Only set constr=False if you remove the intercept from your model. Otherwise, the model is not identifiable and results may be unreliable.
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 350.12 166.34 125.15 318.84 764.26 260.90Higher precision = smoother curve (smaller second differences/curvature). RW2 typically has higher precision values than RW1 due to the different scale of second vs first differences.
Random Effect Estimates
The posterior estimates for each time point are in result.summary_random:
# Smooth temporal trend estimates
trend = result.summary_random['time']
print(trend)
# Plot the smooth curve with uncertainty
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('RW2 Effect')