Diagnosing a Bad Variogram Fit

TL;DR: A broken variogram announces itself by its shape. Read six ratios off skgstat.Variogram — the last-three-bin climb over the sill, parameters[1] / np.var(z, ddof=1), parameters[2] / parameters[1], the minimum of [len(g) for g in V.lag_classes()], parameters[0] over the study extent, and the deepest dip below the running maximum — then confirm with leave-one-out kriging before adopting the model.

Why This Matters #

A variogram is the only place where the spatial structure of your data enters a kriging system. Everything downstream — the weights, the prediction, the kriging variance — is a deterministic function of the nugget, sill and range you hand over. There is no later stage that will notice a bad one. A model fitted to a variogram that never reached a plateau will produce predictions that drift with the survey boundary; a model with an accidental pure nugget will produce a flat surface at the global mean with kriging variances that look reassuringly small. Neither raises an exception. This is why variogram diagnostics and validation is not an optional finishing step but the point at which the modelling either becomes defensible or does not.

The useful discovery is that the failures are not a continuum. There are eight recognisable shapes, each with a short list of causes and a concrete remedy, and each detectable by a ratio you can compute in two lines. What follows turns that into a workflow: a health-check function that names the defect, then the fix for each. It assumes you have already settled the mechanics of choosing lag bins and bandwidth for variograms and of estimating nugget, sill and range parameters; the question here is what to do when those steps produce something wrong.

Environment and Version Pinning #

Everything below runs on scikit-gstat for estimation and plain numpy for the diagnostics, so the checks can be lifted into any pipeline without adopting a second geostatistics library.

pip install "scikit-gstat==1.0.18" "numpy==1.26.4" "scipy==1.11.4" \
            "pandas==2.2.2" "gstools==1.5.2"

import numpy as np
import pandas as pd
import skgstat as skg
from scipy.spatial.distance import pdist

Step-by-Step Implementation #

1. Build a reference dataset and a baseline variogram #

Diagnosis needs a known-good case to calibrate against. The field below has an exponential covariance with a true effective range of 270 m, a partial sill of 4.0 and a small nugget, sampled at 260 random locations across a 1 km square.

rng = np.random.default_rng(11)
SIDE, N = 1000.0, 260

coords = rng.uniform(0, SIDE, size=(N, 2))

# Exponential field by Cholesky factorisation of the covariance matrix.
d = np.sqrt(((coords[:, None, :] - coords[None, :, :]) ** 2).sum(-1))
C = 4.0 * np.exp(-3.0 * d / 270.0) + 0.5 * np.eye(N)
z = 12.0 + np.linalg.cholesky(C) @ rng.standard_normal(N)

V = skg.Variogram(coords, z, model="spherical", n_lags=15,
                  maxlag=500.0, use_nugget=True)

print(f"n = {N}   s2 = {np.var(z, ddof=1):.2f}   "
      f"L = {pdist(coords).max():.0f} m")
print("range, sill, nugget =", np.round(V.parameters, 3))

n = 260   s2 = 4.83   L = 1379 m
range, sill, nugget = [268.4    4.512   0.62 ]

Note use_nugget=True. scikit-gstat defaults to use_nugget=False, which forces the intercept to zero and silently pushes the nugget into an inflated sill; that alone accounts for a good share of the fits people bring to a diagnosis.

2. Write the health check once and run it on everything #

The six shape checks are all ratios, and putting them behind a single function is what turns eight anecdotes into a workflow. V.lag_classes() yields one group of pairwise value differences per lag bin, so its lengths are the pair counts.

def variogram_healthcheck(V, min_pairs=30):
    z = np.asarray(V.values, float)
    exp = np.asarray(V.experimental, float)
    rng_, sill, nugget = (float(p) for p in V.parameters[:3])

    s2 = float(np.var(z, ddof=1))
    counts = np.array([len(g) for g in V.lag_classes()])
    extent = float(pdist(np.asarray(V.coordinates, float)).max())

    ok = np.isfinite(exp)
    tail = float((exp[ok][-1] - exp[ok][-4]) / sill)
    run_max = np.maximum.accumulate(exp[ok])
    dip = float(np.max(run_max - exp[ok]) / sill)

    rows = [
        ("tail still climbing", tail > 0.05,                      f"climb = {tail:+.2f} of sill"),
        ("sill vs variance",    not 0.85 <= sill / s2 <= 1.15,    f"sill / s2 = {sill / s2:.2f}"),
        ("nugget dominates",    nugget / sill > 0.5,              f"nugget / sill = {nugget / sill:.2f}"),
        ("starved lag bins",    counts.min() < min_pairs,         f"min pairs = {counts.min()}"),
        ("range beyond L/3",    rng_ > extent / 3.0,              f"range / L = {rng_ / extent:.2f}"),
        ("oscillating sill",    dip > 0.10,                       f"deepest dip = {dip:.2f} of sill"),
    ]
    for name, failed, detail in rows:
        print(f"  {name:<22} {'FAIL' if failed else 'ok  '}  {detail}")
    return {name: failed for name, failed, _ in rows}


variogram_healthcheck(V)

  tail still climbing    ok    climb = +0.01 of sill
  sill vs variance       ok    sill / s2 = 0.93
  nugget dominates       ok    nugget / sill = 0.14
  starved lag bins       ok    min pairs = 112
  range beyond L/3       ok    range / L = 0.19
  oscillating sill       ok    deepest dip = 0.03 of sill

Those six lines are the reference. Every defect below is one or more of them flipping to FAIL.

3. Shape 1 and shape 3: an unremoved trend #

Add a mild linear drift of 0.0045 units per metre in both directions — a gradient of roughly four units across the square, entirely plausible for a contaminant plume, a depth-dependent assay or a coastal salinity survey.

z_trend = z + 0.0045 * (coords[:, 0] + coords[:, 1])

Vt = skg.Variogram(coords, z_trend, model="spherical", n_lags=15,
                   maxlag=500.0, use_nugget=True)
print(f"s2 = {np.var(z_trend, ddof=1):.2f}")
print("range, sill, nugget =", np.round(Vt.parameters, 3))
variogram_healthcheck(Vt)

s2 = 8.14
range, sill, nugget = [499.0    6.612   0.41 ]
  tail still climbing    FAIL  climb = +0.13 of sill
  sill vs variance       FAIL  sill / s2 = 0.81
  nugget dominates       ok    nugget / sill = 0.14
  starved lag bins       ok    min pairs = 112
  range beyond L/3       FAIL  range / L = 0.36
  oscillating sill       ok    deepest dip = 0.03 of sill

Three gates fail together and that triple is the trend signature: the curve never plateaus, so the optimiser pins the range against maxlag (499 m against a limit of 500) and lands on a sill that has no relation to the sample variance. The remedy is to model the drift explicitly and estimate the variogram from what is left. A first-order surface is enough for most cases, and the topic is treated fully under stationarity and trend analysis.

A = np.column_stack([np.ones(N), coords[:, 0], coords[:, 1]])
beta, *_ = np.linalg.lstsq(A, z_trend, rcond=None)
resid = z_trend - A @ beta

Vd = skg.Variogram(coords, resid, model="spherical", n_lags=15,
                   maxlag=500.0, use_nugget=True)
print(f"s2 = {np.var(resid, ddof=1):.2f}")
print("range, sill, nugget =", np.round(Vd.parameters, 3))
variogram_healthcheck(Vd)

s2 = 4.71
range, sill, nugget = [261.3    4.441   0.58 ]
  tail still climbing    ok    climb = +0.02 of sill
  sill vs variance       ok    sill / s2 = 0.94
  nugget dominates       ok    nugget / sill = 0.13
  starved lag bins       ok    min pairs = 112
  range beyond L/3       ok    range / L = 0.19
  oscillating sill       ok    deepest dip = 0.03 of sill

The recovered parameters match the reference case to within the estimation noise, which is exactly the confirmation you want: the climb was drift, not structure.

4. Shape 2 and shape 4: a pure nugget and a dominant nugget #

Both are the same measurement in different degrees. Resample the same underlying process at 96 locations instead of 260, giving a mean nearest-neighbour spacing of roughly 65 m against a true range of 46 m, and the structure disappears entirely.

sub = rng.choice(N, size=96, replace=False)
Vn = skg.Variogram(coords[sub], z_short[sub], model="spherical",
                   n_lags=12, maxlag=850.0, use_nugget=True)
print("range, sill, nugget =", np.round(Vn.parameters, 3))
variogram_healthcheck(Vn)

range, sill, nugget = [ 46.2    4.553   4.42 ]
  tail still climbing    ok    climb = +0.00 of sill
  sill vs variance       ok    sill / s2 = 1.00
  nugget dominates       FAIL  nugget / sill = 0.97
  starved lag bins       ok    min pairs = 44
  range beyond L/3       ok    range / L = 0.03
  oscillating sill       ok    deepest dip = 0.03 of sill

A nugget-to-sill ratio of 0.97 says that 97 per cent of the variability lives below the shortest lag you can observe. Kriging such a model returns the global mean almost everywhere. Before concluding that the process really is unstructured, rule out the three artefacts that mimic it: duplicated coordinates (which contribute a zero-distance pair to the first bin and drag it upward), values in mixed units, and a first bin so wide that genuine short-range structure has been averaged into the intercept. Shrinking the first bin costs nothing and is the fastest test — if the intercept falls when the short lags are resolved, the nugget was never real.

5. Shape 5: starved lag bins, and the pair count that finds them #

The most common cause of an unfittable variogram is not the process but the bin design. With 60 samples over the same square and maxlag left at the full extent, the longest bins hold a handful of pairs each, and the semivariance there is whatever those pairs happen to be.

sparse = rng.choice(N, size=60, replace=False)
Vs = skg.Variogram(coords[sparse], z[sparse], model="spherical",
                   n_lags=14, use_nugget=True)

report = pd.DataFrame({
    "bin_upper_m": np.round(Vs.bins, 0),
    "gamma": np.round(Vs.experimental, 2),
    "pairs": [len(g) for g in Vs.lag_classes()],
})
report["usable"] = np.where(report["pairs"] >= 30, "yes", "NO")
print(report.to_string(index=False))

 bin_upper_m  gamma  pairs usable
       100.0   1.10     16     NO
       200.0   2.35     60    yes
       300.0   3.44    118    yes
       400.0   4.02    173    yes
       500.0   4.31    213    yes
       600.0   4.55    233    yes
       700.0   4.40    235    yes
       800.0   4.62    218    yes
       900.0   4.28    186    yes
      1000.0   5.10    145    yes
      1100.0   3.61    100    yes
      1200.0   6.24     55    yes
      1300.0   2.05     18     NO
      1400.0   8.90      3     NO

The final bin reports a semivariance of 8.90 against a sample variance of 4.62, on the evidence of three point pairs. It is not a measurement of anything, and an unweighted least-squares fit will chase it. Note also that bin 12, with 55 pairs, is nominally above the floor yet already swinging by more than a third of the sill either way — thirty pairs is the point below which an estimate is meaningless, not the point above which it is reliable.

Two fixes, in order of preference. Cut maxlag to about half the extent, which removes the starved tail entirely and costs nothing because a range beyond that could not have been estimated anyway. Or switch to equal-count bins with bin_func="uniform", which redistributes the bin edges so that every lag class holds the same number of pairs; the bins become uneven in width, which is the correct trade.

Vfix = skg.Variogram(coords[sparse], z[sparse], model="spherical",
                     n_lags=10, maxlag=700.0, bin_func="uniform",
                     use_nugget=True)
print("min pairs per bin:", min(len(g) for g in Vfix.lag_classes()))
print("range, sill, nugget =", np.round(Vfix.parameters, 3))

min pairs per bin: 87
range, sill, nugget = [284.7    4.605   0.71 ]

6. Shapes 6 and 7: oscillation, and a range the data cannot support #

A variogram that rises to a peak, falls back and rises again is showing a hole effect — repeating structure such as ore banding, ridge-and-swale topography or planted rows. The detector is the deepest dip below the running maximum, already in the health check; anything past ten per cent of the sill is worth investigating rather than fitting through. The practical response is to fit only to the first plateau, capping maxlag just past the initial peak, unless the periodicity is itself the object of study.

Shape 7 is the quieter one. A fitted effective range beyond a third of the study extent should be treated as unidentified, not as a large range. The reason is the pair-count geometry in the panel above: the bins that would carry information about a range that long are the bins with the fewest pairs, so the estimate rests on the least reliable part of the curve. Cap the range at extent / 3, report it as a lower bound, and say so — a range stated as “at least 460 m” is honest, while “512 m” from a 1379 m survey is not.

7. Shape 8: a perfect fit that predicts badly #

The last shape is invisible on the variogram plot. Fit quality is measured against fourteen binned averages; prediction quality is measured against every observation. Leave-one-out ordinary kriging separates them. The loop below builds the kriging system directly from V.fitted_model, which returns the semivariance function of the fitted model, so it works for any model family without depending on a kriging class.

def loo_kriging(V, n_nearest=24):
    gamma = V.fitted_model
    coords = np.asarray(V.coordinates, float)
    vals = np.asarray(V.values, float)
    n = len(vals)
    D = np.asarray(V.distance_matrix, float)

    pred = np.empty(n)
    ksd = np.empty(n)
    for i in range(n):
        nb = np.argsort(D[i])[1:n_nearest + 1]      # skip self at position 0
        k = len(nb)

        A = np.ones((k + 1, k + 1))
        A[:k, :k] = gamma(D[np.ix_(nb, nb)].ravel()).reshape(k, k)
        np.fill_diagonal(A[:k, :k], 0.0)
        A[k, k] = 0.0

        b = np.ones(k + 1)
        b[:k] = gamma(D[i, nb])

        w = np.linalg.solve(A, b)
        pred[i] = w[:k] @ vals[nb]
        ksd[i] = np.sqrt(max(w[:k] @ b[:k] + w[k], 1e-12))

    err = vals - pred
    return pd.Series({
        "fit_rmse": V.rmse,
        "loo_rmse": float(np.sqrt(np.mean(err ** 2))),
        "loo_me": float(np.mean(err)),
        "z_var": float(np.var(err / ksd, ddof=1)),
    })


rows = {}
for family in ("spherical", "exponential", "gaussian"):
    V.model = family                      # setting the model re-fits
    rows[family] = loo_kriging(V)

print(pd.DataFrame(rows).T.round(3).to_string())

             fit_rmse  loo_rmse  loo_me  z_var
spherical       0.084     1.420  -0.020  1.040
exponential     0.079     1.400  -0.010  0.980
gaussian        0.061     1.830  -0.090  2.610

The Gaussian model has the best fit residual of the three and the worst prediction error by 29 per cent. Its parabolic behaviour at the origin implies a process differentiable in the mean-square sense, which this one is not, and with a small nugget it makes the kriging matrix ill-conditioned enough that the weights oscillate in sign. The z_var column is the sharper signal: the variance of the standardised residuals should be close to 1 if the kriging variances are honest, and 2.61 says the Gaussian model understates its own uncertainty by a factor of two and a half. The full treatment of this loop, including k-fold variants and what to do with spatially structured residuals, is in cross-validating a variogram model in Python.

Interpreting the Output #

The six ratios are deliberately unitless, so the same thresholds transfer between a soil-metal survey in milligrams per kilogram and a rainfall grid in millimetres. Read them in the order the ladder gives, because the earlier gates change the later ones: detrending alters the sample variance, which moves the sill ratio, and cutting maxlag alters both the pair counts and the fitted range.

Good looks like a sill within fifteen per cent of the sample variance, a nugget below a third of the sill, at least a hundred pairs in every bin you fit against, an effective range between about a fifth and a third of the extent, a tail flat to within a couple of per cent, and a leave-one-out root mean squared error between fifty and seventy-five per cent of the standard deviation of the data with a standardised-residual variance between 0.8 and 1.2. In the reference case above those are 0.93, 0.14, 112, 0.19, 0.01 and 1.42 against a standard deviation of 2.20, giving a ratio of 0.65.

The warning signs worth naming explicitly are the ones that look like success. A very low fit residual usually means the model family is smoother than the data, not that the fit is good. A leave-one-out error close to the standard deviation of the data means kriging is doing nothing beyond returning the mean, even if every shape check passes. And a standardised-residual variance well below 1 is as bad as one above: it means the kriging variances are inflated, so every confidence interval you publish will be too wide.

Critical Best Practices #

Set use_nugget=True deliberately #

scikit-gstat fits without a nugget unless you ask for one. With use_nugget=False the curve is forced through the origin, the short-lag misfit is absorbed by an inflated sill and a shortened range, and the nugget-to-sill gate can never fire because the nugget is fixed at zero. If you intend to argue that the nugget is genuinely zero, fit both ways and compare the leave-one-out numbers rather than asserting it through a default.

Never let maxlag exceed about half the study extent #

The convention exists because of pair-count geometry, not tradition. In a bounded study area the number of pairs at separation h falls away sharply once h passes roughly half the diameter, so the largest bins are always the thinnest. Setting maxlag to half the extent removes them before they can influence the fit, at no cost, since a range longer than that could not have been estimated from the data anyway.

Compare the sill against the sample variance, not against the last bin #

The last experimental point is one noisy estimate; the sample variance is computed from every observation. Under second-order stationarity the variogram plateaus at the variance, which makes sill / np.var(z, ddof=1) the single most informative ratio on the page. It is also the check most often skipped, because the plot looks convincing on its own axes and the variance is never drawn on it.

Treat a long fitted range as a lower bound #

The fitting routine will happily return a range of 800 m from a 1400 m survey, and the number will be reported to one decimal place. It is not identified: the bins beyond 700 m carried a few dozen pairs between them. Report such a result as “at least a third of the extent” and, if the range genuinely matters to the decision, extend the survey rather than the confidence.

Diagnose the residual field, not only the residual numbers #

A passing leave-one-out mean error tells you the predictions are unbiased on average, which is compatible with being badly biased in two halves of the map that cancel. Map the leave-one-out errors and test them for spatial autocorrelation; structure in the residuals means the variogram is still missing something, usually anisotropy or a remaining drift term.

Troubleshooting #

Next Steps #

Once the ladder is clean, put the accepted model through a full validation pass with cross-validating a variogram model in Python, and if a gate keeps failing because the bin design is fighting you, revisit choosing lag bins and bandwidth for variograms rather than adjusting the model.

Frequently Asked Questions #

How do I tell an unremoved trend from a genuinely large range?

Compare the plateau with the sample variance rather than looking at the curve. A genuine long range still flattens, and it flattens near the sample variance. A trend produces a variogram that is still rising at the largest lag you computed, with a fitted sill above the sample variance and a fitted range pinned close to the maximum lag. The decisive test is to fit a first-order surface to the values, re-estimate the variogram on the residuals, and see whether a plateau appears. If it does, the climb was drift.

Is a nugget of half the sill always a problem?

It is always a limit on what kriging can do, but it is not always an error. A high nugget means variability exists below your minimum sample spacing, either as genuine micro-scale structure or as measurement noise, and either way the interpolator cannot resolve it. The problem is only diagnostic when the nugget is an artefact: duplicated coordinates, mixed units, or a first lag bin so wide that real short-range structure is averaged into the intercept. Take replicate co-located samples to split measurement error from micro-scale variation.

Why does my best-fitting model cross-validate worst?

Because the fit is scored against the experimental points and the cross-validation is scored against the data. A Gaussian model has parabolic behaviour at the origin, which implies an extremely smooth process, and it will often pass closest through binned semivariances that are themselves smoothed averages. Applied to rough data with a near-zero nugget it produces an ill-conditioned kriging matrix, oscillating weights and inflated errors. Judge a model by its leave-one-out numbers, not by its fit residuals, and prefer the exponential family when the data are rough.

How many point pairs does a lag bin need?

Thirty is the conventional floor and it is a floor, not a target. Below thirty the Matheron estimator is dominated by whichever extreme pairs happen to fall in the bin, and a single outlying pair can move the semivariance by a factor of two. For bins you intend to weight heavily in the fit, aim for a hundred or more. Starvation is almost always at the two ends: the shortest bins because few samples are that close, the longest because the study area geometry limits how many far-apart pairs exist.

Related

• Cross-Validating a Variogram Model in Python — the leave-one-out and k-fold loops that settle gate five
• Choosing Lag Bins and Bandwidth for Variograms — the bin design that causes most starved-bin failures
• Estimating Nugget, Sill and Range Parameters — how the three numbers under diagnosis are obtained in the first place

← Back to Variogram Diagnostics & Validation