MLE Cross-Experiment Comparison Guide

Overview

Standard validation metrics (RMSE, bias, correlation) compare each experiment against observations in isolation. When you have three or more experiments that share the same observation dataset — sensitivity runs, ensemble members, or physics variants — a scalar metric gives no information about the joint structure of the errors.

MLE comparison in reduced-rank (SVD) space addresses this. The joint residual matrix across all experiments is decomposed with SVD, and each experiment is ranked by its Gaussian log-likelihood in the leading modes. This is the same “optimal fingerprinting” framework used in climate attribution (Allen & Tett 1999): the leading SVD modes capture the dominant joint variance pattern, and the log-likelihood in that low-dimensional space is a principled, multi-dimensional distance from perfect.

The result is a ranking table where the best experiment has Δlog L = 0 and all others have Δlog L < 0. A large negative Δlog L means that experiment’s residual pattern is substantially different from the best, in the directions that matter most across the ensemble.

Why SVD?

The residual matrix R has shape n_obs × n_exp (often 100 000 × 8). A naive multivariate Gaussian log-likelihood in the full n_obs-dimensional space is intractable for two reasons:

1. The covariance matrix cannot be inverted. With only n_exp columns, the sample covariance of R has rank at most n_exp, far smaller than n_obs. The inverse does not exist and the Gaussian likelihood is undefined.

2. The column space of R is at most n_exp-dimensional. No matter how many observations there are, experiments can only differ along at most n_exp independent directions. Anything outside that subspace is shared error and contributes nothing to the ranking.

SVD solves both problems by working exclusively in the low-dimensional subspace spanned by the experiment columns. The decomposition R_std = U diag(s) Vᵀ reveals that subspace via the left singular vectors U. Keeping the r leading modes (by the variance threshold) retains the directions of maximum joint variance — the directions where experiments most clearly differ — and discards the noise-dominated tail.

Within this r-dimensional space the Gaussian log-likelihood is well-posed: there are at most n_exp − 1 non-trivial modes, r is small, and the remaining variance (the discarded tail) is absorbed into the σ² estimate. The result is a principled, noise-regularised distance from perfect that is comparable across experiments regardless of how many observations exist.

A secondary benefit is numerical: the SVD of an n_obs × n_exp matrix costs O(n_obs × n_exp²), which is fast even for large observation datasets because n_exp is always small (see Computational cost).

When to use this

• You have three or more experiments with different physics, forcing, or grid settings.
• All experiments were evaluated against the same observation dataset (e.g. Argo floats or GLODAP profiles).
• The experiments cover overlapping time periods and geographic domains.

It is less useful when:

• You only have two experiments — the result collapses to a weighted RMSE comparison and provides no additional information.
• Experiments cover non-overlapping years — the common observation footprint will be empty or too small.

How it works

experiment DataFrames (time, lat, lon, depth, value, model_value, variable)
         │
         ▼
1. Common footprint
   Intersect observation keys present in ALL experiments with non-NaN obs and model values.

         │
         ▼
2. Residual matrix R  [n_obs × n_exp]
   R[i, k] = model_value[i,k] − obs_value[i]

         │
         ▼
3. Standardise per variable
   Divide each variable's rows by the obs standard deviation so TEMP, PSAL, and DOXY
   are commensurable.

         │
         ▼
4. SVD  →  R_std ≈ U · diag(s) · Vt
   Select leading r modes by cumulative variance threshold (default 90%).

         │
         ▼
5. Project each experiment: r̃_k = U_r^T · R_std[:, k]   shape (r,)

         │
         ▼
6. Log-likelihood: log L_k = −½ σ⁻² ‖r̃_k‖²
   Rank by Δlog L = log L_k − max(log L)

         │
         ▼
7. Additive breakdown (diagnostics only)
   For each group g (variable or obs type), rows are disjoint so:
       r̃_k = Σ_g  U_r[g_mask]ᵀ R_std[g_mask, k]   (exact)
   Per-group log-likelihood contribution:
       log L_{k,g} = −½ σ⁻² r̃_{k,g}ᵀ r̃_k
   Σ_g log L_{k,g} = log L_k  (exact, no cross-term correction needed)

Prerequisites

Each experiment must first produce a residuals parquet file. The path and filename depend on the analysis type:

Run each analysis once per experiment, changing experiment: in the config. All experiments must share the same area: so the observation footprint is comparable.

CLI usage

Step 1: Generate residual files

# Profile observations (save_residuals: true is the default)
argo-profiles --config config/baseline.yaml
argo-profiles --config config/sensitivity1.yaml
argo-profiles --config config/sensitivity2.yaml

# Gridded surface obs (add save_residuals: true + spatial_thinning to the YAML)
gridded-2d-validation --config config/baseline_2d.yaml
gridded-2d-validation --config config/sensitivity1_2d.yaml

# Gridded 3D obs (same)
gridded-3d-validation --config config/baseline_3d.yaml
gridded-3d-validation --config config/sensitivity1_3d.yaml

Step 2: Run MLE comparison

# Auto-discover all experiments with Argo residuals for this area
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types argo --variables TEMP PSAL DOXY --diagnostics

# Exclude a specific experiment from auto-discovery
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types argo --exclude-experiments OldRun TestRun

# Explicitly name experiments (disables auto-discovery)
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --experiments Baseline Sensitivity1 Sensitivity2 \
    --obs-types argo --variables TEMP PSAL DOXY

# Combine Argo (deep water) + fixed stations (coastal) for fuller coverage
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types argo fixed_platform --variables TEMP PSAL

# Gridded surface obs (e.g. satellite SST, bottom-T reanalysis)
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types surface --variables TEMP

# Combine gridded surface obs with Argo for joint surface + profile ranking
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types surface argo --variables TEMP PSAL

# Gridded 3D obs (e.g. WOA23 climatology)
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types gridded3d --variables TEMP PSAL

# Log-transform phytoplankton and oxygen before ranking
mle-comparison --analyses-dir $OCEANICU_ANALYSES_FOLDER --area NS \
    --obs-types argo --variables TEMP PSAL CHLA DOXY \
    --transforms CHLA:log DOXY:log

When --experiments is omitted the script scans analyses/areas/<area>/validations/ for every experiment directory that contains <obs_type>_residuals.parquet for all requested obs types. The mle_comparison output directory is always skipped.

Residuals from all --obs-types are concatenated per experiment before the common footprint is built, so each obs source contributes its locations.

Key options:

Variable transforms

For biogeochemical variables whose distributions are right-skewed (phytoplankton chlorophyll, nutrients, turbidity), a Gaussian residual assumption is not appropriate in linear space. Applying a log-transform before the MLE makes the residuals approximately Gaussian:

Rows where the raw value or model value is ≤ 0 are dropped before the common footprint is built (a warning is logged for each experiment).

Using a config file

# config/mle_comparison.yaml
area: NS
analyses_dir: /data/local/analyses   # a real path — ${VAR} is NOT expanded here
# experiments: [Baseline, Sensitivity1, Sensitivity2]  # omit to auto-discover
exclude_experiments: [OldRun]   # dropped from auto-discovery (ignored if experiments is set)
obs_types: [argo, fixed_platform]   # list for multi-source, or a single string
variables: [TEMP, PSAL, CHLA]
transforms:
  CHLA: log       # ln(model/obs) residuals for chlorophyll
  DOXY: log       # ln(model/obs) residuals for dissolved oxygen
rank: 3
diagnostics: true

mle-comparison --config config/mle_comparison.yaml

Command-line flags override config file values.

Output

Results are written to analyses/areas/<area>/validations/mle_comparison/:

mle_ranking.csv columns:

Why gridded validation scores may diverge from MLE rankings

It is common to see a large improvement in a gridded 3D or surface RMSE score for one experiment while the MLE ranking shows little or no separation. This is not a bug; it reflects a genuine difference in what the two metrics measure.

Gridded validation (e.g. gridded-3d-validation) compares the model to a gridded reference product — a WOA climatological mean, a CMEMS reanalysis, or a satellite analysis. These products are spatially dense (every grid cell has a value) but heavily smoothed and, in the case of reanalyses, data-assimilated. An experiment that is initialised from or constrained by the same data family will naturally score better on that metric. The score also measures bulk domain-mean skill averaged over every grid cell, which can be dominated by a few large regions or seasons.

MLE ranking uses raw in-situ residuals from ARGO floats or ICES profiles — sparse, independent point measurements at actual ocean interior locations with no data assimilation applied. It asks: does this experiment produce residuals that are collectively closer to Gaussian noise at the places where observations exist? This is a stricter and more independent test.

The divergence is therefore informative:

Use both metrics together. Neither supersedes the other; they answer different questions about model skill.

Interpreting results

Only delta_log_L matters, not the absolute log-likelihood value. delta_log_L = 0 identifies the best experiment; more negative values indicate larger residuals in the joint SVD space.

With two experiments the result is equivalent to a weighted RMSE comparison. The method adds value from three experiments upward.

AIC (−2 log L + 2) can be used in a model-selection sense: prefer the experiment with the lowest AIC. The + 2 term penalises the single free parameter σ². For the purposes of cross-experiment ranking, AIC and Δlog L give the same ordering.

Scree plot

Shows the singular values and their cumulative explained variance. If a single mode dominates, the ranking is essentially one-dimensional and driven by one error pattern. If several modes have similar singular values, the comparison is genuinely multi-dimensional and the mode-loading plots for each mode are all relevant.

Mode-loading plot

Each SVD mode has a corresponding left singular vector U of length n_obs. The mode-loading plot maps the top 30 entries by magnitude onto geographic space, where:

• Circle size is proportional to the absolute value of the loading — larger circles are the observation locations that most strongly define this mode and carry the most weight in separating experiments.
• Circle colour is the signed loading value. Locations with the same colour move together in this mode: when one experiment has a larger residual at a dark-red location, it tends to have larger residuals at all other dark-red locations too. Locations with opposite colours (red vs. blue) are anti-correlated in this mode.
• Labels show the variable (e.g. PSAL, TEMP) at that location.

What to look for:

• All circles the same colour → mode 0 captures a nearly uniform bias that all experiments share. The ranking reflects which experiment has the smallest overall bias in that direction.
• A mix of red and blue → the mode captures a spatial contrast (e.g. north–south or coastal–offshore gradient). Experiments that do well in one region but poorly in another will be separated here.
• Circles concentrated in one sub-region → that geographic area is driving the ranking. Check whether the model resolution or forcing is different there.

Example: if the mode-0 map shows all-negative loadings of similar magnitude spread across the domain (as in a uniform salinity bias), mode 0 is telling you that the experiments differ mainly in their domain-wide mean salinity error, not in their spatial pattern. An experiment with a smaller mean salinity bias will have a less-negative log-likelihood contribution in this mode.

Per-mode contributions bar chart

Shows each experiment’s log-likelihood contribution broken down by SVD mode. One group of bars per mode (left to right), one bar per experiment per group.

• Bar height (all negative) is −½ σ⁻² r̃² for that experiment in that mode: a shorter bar means the residual projected onto this mode is smaller, i.e. better.
• When only one mode is shown (rank = 1), the total ranking collapses to a single number. All experiments are compared on the same dominant error pattern and their relative bar heights directly mirror the delta_log_L ranking table.
• When multiple modes are shown, you can see whether an experiment that wins overall does so because it is good across all modes, or because it excels in one mode while being mediocre in others.

Example: if seven experiments all show similar bar heights in a single-mode result (bars clustered between −0.3 and −0.5), the experiments are not strongly separated. A large spread (one bar near 0, another near −1) indicates a clear winner. Bars of almost identical height mean the observation dataset does not have enough discriminating power to separate those two experiments — consider adding more obs types or variables.

Contributions by variable / obs type

These charts decompose the total log-likelihood into additive per-variable (or per-obs-type) contributions. The decomposition is exact: because the observation rows are partitioned by group, the group projected residuals sum to the total:

r̃_k = Σ_g  U_r[g_mask]ᵀ R_std[g_mask, k]      (exact, not approximate)

The additive contribution of group g to experiment k’s log-likelihood is then −½ σ⁻² r̃_{k,g}ᵀ r̃_k, and summing over groups recovers the total log L_k = −½ σ⁻² ‖r̃_k‖² exactly. The bars therefore stack to the total log-likelihood bar for each experiment.

• Bars shorter (less negative) for a group → that variable or obs source is where the experiment performs best relative to the others.
• One group dominates all experiments equally → the ranking is insensitive to that group; removing it would not change the order.
• One group separates the experiments → that variable or obs source carries the discriminating power. Cross-check with the mode-loading map to see which geographic region is responsible.

Rank selection

The default automatic selection combines two criteria:

1. Variance threshold: smallest r where cumulative explained variance ≥ variance_threshold (default 0.90).
2. Elbow: the r at the largest second-difference drop in the singular value spectrum.

The chosen rank is min(threshold_r, elbow_r), clamped to [1, n_modes − 1].

Override with --rank when you have domain knowledge about how many independent error patterns to expect. Too low a rank discards real signal; too high a rank projects onto noise and compresses the ranking. The scree plot is the diagnostic to consult.

Programmatic use

from lib.mle_comparison import compare_experiments_mle

# experiment_dfs: dict of experiment_name -> DataFrame
result = compare_experiments_mle(
    experiment_dfs,
    variables=['TEMP', 'PSAL', 'DOXY'],
    variance_threshold=0.90,
    return_diagnostics=True,
)
ranking_df, diag = result

print(ranking_df)
# Diagnostics dict keys:
#   singular_values, explained_variance, chosen_rank,
#   U_r, Vt_r, r_tilde, sigma2, meta_df, scalings

Each DataFrame in experiment_dfs must have columns: time, lat, lon, depth, value, model_value, variable and optionally profile_id.

Scope and the independence assumption

The log-likelihood formula assumes independent observations:

log L_k = −½ σ⁻² ‖r̃_k‖²

This is naturally satisfied for sparse in-situ data (Argo, GLODAP, WOD, cruise CTD) where floats and casts are separated by tens to hundreds of kilometres.

Gridded data — use with thinning

gridded-2d-validation and gridded-3d-validation can also write residuals parquets (see those guides). Dense gridded data violates the independence assumption: a 0.1° North Sea grid might have 50,000 points but a spatial correlation length of ~50 km, giving perhaps a few hundred effective degrees of freedom. Without thinning, Δlog L values will look larger and more significant than they truly are.

The spatial_thinning config in those scripts mitigates this:

output:
  save_residuals: true
  spatial_thinning:
    stride: 5         # keep every 5th grid point in lat and lon
    z_stride: 2       # (3D only) keep every 2nd depth level
    max_rows: 100000  # hard cap

Practical guidance:

• Start with stride: 5 and max_rows: 100000. Re-run mle-comparison with stride: 3 and stride: 8; if the relative ranking is stable, the thinning level is sufficient.
• For 3D gridded obs, also use z_stride: 2 or z_stride: 3 to reduce depth-level correlation.
• Combining thinned gridded obs with profile obs (--obs-types surface argo) is often more informative than either alone: gridded obs provide dense spatial coverage of the surface; profiles provide depth structure at sparse locations.
• The SVD step is the same regardless of obs type — the leading modes capture the dominant joint-variance pattern across all supplied obs.

For the most defensible results where spatial covariance structure matters, average to regional boxes first and apply MLE to the box means, where the independence assumption is much more defensible.

Computational cost

The method is not computationally intensive. The matrix being decomposed is n_obs × n_exp, and n_exp is always small (typically 3–10).

The SVD is the mathematically dominant step, but because n_exp is tiny scipy uses the thin SVD and the cost scales linearly with n_obs. A 100,000 × 5 matrix takes well under a second.

The practical bottleneck is reading and joining the parquet files, not the maths. The residual matrix R lives entirely in RAM: 1 million obs × 10 experiments × 8 bytes = 80 MB, which is negligible for realistic datasets.

For gridded data with max_rows: 100000 and stride: 5, the combined parquet stays well under 100 MB and the SVD completes in seconds. Without thinning (full-resolution gridded output) n_obs could reach tens of millions and the in-memory matrix would become expensive; always apply spatial_thinning (see Scope section above).

Limitations

• Common obs footprint required: all experiments must cover the same stations and years. Experiments with non-overlapping periods cannot be compared.
• Diagonal covariance: observation errors are treated as independent and equal within each variable. Spatial correlation between obs is ignored. This is a reasonable approximation for scattered profile data and for thinned gridded data; it becomes increasingly problematic as thinning is relaxed (see Scope above).
• Low power with few experiments: with fewer than four experiments the SVD modes are poorly constrained and rankings may be unreliable.
• Pooled σ²: the noise variance is estimated across all experiments. If one experiment is catastrophically wrong, it inflates σ² and compresses the Δlog L differences among the others.