Regression Kriging vs Universal Kriging

TL;DR: With a linear drift and the same variogram, regression kriging and universal kriging are the same estimator. gstools.krige.ExtDrift(model, cond_pos, cond_val, ext_drift=elev) reproduces a GLS trend fit plus simple kriging of the residuals to machine precision. The differences that matter are that ordinary least squares biases the residual variogram low, and that naive regression kriging omits the trend-uncertainty term from its variance.

Why This Matters #

The two methods are usually presented as rivals, and a good deal of ink has been spent on which predicts better. The honest answer is that when the trend is linear in its coefficients and both use the same covariance model, they predict identically — not similarly, identically, to twelve decimal places. The result is a corollary of the best linear unbiased predictor: universal kriging is generalised least squares on the trend plus simple kriging on what is left, rearranged into a single linear system. So the interesting question is never “which is more accurate” but “which failure modes am I buying”.

That reframing changes what you check. The regression kriging route wins whenever the trend needs a model class kriging cannot express — a random forest, a boosted tree ensemble, a generalised additive model — and it lets you look at the residual variogram directly, which is the only honest way to check that the trend has left something stationary behind. The ordinary and universal kriging route wins whenever the drift really is linear, because it gives you a variance that already accounts for not knowing the trend coefficients, and there is no bookkeeping to get wrong. Both depend on the same prior judgement about what is trend and what is residual, which is the subject of stationarity and trend analysis.

The Algebra, Stated Once #

Write the model as $Z = F\beta + e$, where $F$ is the $n \times p$ matrix of drift terms evaluated at the sample locations (a column of ones and a column of elevations, say), $\beta$ the unknown coefficients, and $e$ a zero-mean residual with covariance matrix $C$. Let $c_0$ hold the covariances between the target location $s_0$ and the $n$ samples, $f_0$ the drift terms at $s_0$, and $\sigma_0^2$ the total variance at a point.

The generalised least squares estimate of the trend is

$$\hat{\beta} = (F^{\mathsf{T}} C^{-1} F)^{-1} F^{\mathsf{T}} C^{-1} z ,$$

and the best linear unbiased predictor is

$$\hat{Z}(s_0) = f_0^{\mathsf{T}} \hat{\beta} + c_0^{\mathsf{T}} C^{-1} (z - F\hat{\beta}) .$$

Read left to right, that second equation is regression kriging: a fitted trend surface plus a simple-kriging interpolation of the residuals. Rearranged into a single system with Lagrange multipliers enforcing $\sum_i \lambda_i f_k(s_i) = f_k(s_0)$, it is universal kriging. The prediction variance of that predictor is

$$\sigma^2_{\mathrm{UK}}(s_0) = \underbrace{\sigma_0^2 - c_0^{\mathsf{T}} C^{-1} c_0}_{\text{residual kriging term}} + \underbrace{d_0^{\mathsf{T}} (F^{\mathsf{T}} C^{-1} F)^{-1} d_0}_{\text{trend uncertainty}} , \qquad d_0 = f_0 - F^{\mathsf{T}} C^{-1} c_0 .$$

Naive regression kriging reports only the first term. The second is non-negative and vanishes only where the data pin the trend down completely, so the naive variance is always an underestimate — never an overestimate, and never equal except in a limit.

Environment and Version Pinning #

Two kriging libraries appear below because each is clearer for one half of the comparison: gstools for the drift systems, scikit-gstat for empirical variogram fitting. pykrige is used only as an independent check on the universal kriging solution.

pip install "numpy==1.26.4" "scipy==1.13.1" "pandas==2.2.2" \
            "gstools==1.5.2" "scikit-gstat==1.0.18" "pykrige==1.7.2" \
            "scikit-learn==1.5.1"

import numpy as np
import gstools as gs
import skgstat as skg
from scipy.spatial.distance import pdist, squareform
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold

Step-by-Step Implementation #

1. Simulate a survey with a known trend and a known residual structure #

Everything below is checkable only because the truth is known: soil organic carbon rises linearly with elevation, and the departure from that line is an exponential random field with an effective range of 1500 m plus a nugget.

rng = np.random.default_rng(2026)
N, SIDE = 180, 12_000.0

def elevation(x, y):
    """A smooth deterministic DEM, known everywhere on the grid."""
    return 80.0 + 0.014 * x + 0.010 * y + 30.0 * np.sin(x / 2600.0)

xs, ys = rng.uniform(0, SIDE, N), rng.uniform(0, SIDE, N)
elev_s = elevation(xs, ys)

# Correlated residual: gstools len_scale is the e-folding scale,
# so an effective range of 1500 m means len_scale = 500 m.
truth = gs.Exponential(dim=2, var=0.35, len_scale=500.0)
e_corr = gs.SRF(truth, seed=20260807)((xs, ys))

soc = 1.20 + 0.0085 * elev_s + e_corr + rng.normal(0.0, np.sqrt(0.05), N)
print(f"elevation {elev_s.min():.0f}-{elev_s.max():.0f} m, "
      f"SOC {soc.min():.2f}-{soc.max():.2f} %")

elevation 113-330 m, SOC 1.11-5.02 %

2. Fit the trend by OLS and look at the residual variogram #

This is the textbook regression kriging first stage, and it is where the first real difference appears.

coords = np.column_stack([xs, ys])
F = np.column_stack([np.ones(N), elev_s])

beta_ols, *_ = np.linalg.lstsq(F, soc, rcond=None)
r_ols = soc - F @ beta_ols

V = skg.Variogram(coords, r_ols, model="exponential",
                  n_lags=15, maxlag=6000, normalize=False)
rng_o, psill_o, nug_o = V.parameters      # effective range, partial sill, nugget
print(f"OLS   b0={beta_ols[0]:.4f}  b1={beta_ols[1]:.6f}")
print(f"resid variogram: range={rng_o:.0f} m  psill={psill_o:.4f}  "
      f"nugget={nug_o:.4f}  total sill={psill_o + nug_o:.4f}")

OLS   b0=0.9817  b1=0.009214
resid variogram: range=1119 m  psill=0.2712  nugget=0.0578  total sill=0.3290

The truth was a range of 1500 m and a total sill of 0.400. The OLS-residual variogram is 25% short on range and 18% low on sill. This is not sampling noise: least squares that ignores the correlation between residuals absorbs part of the long-range correlated variation into the fitted trend, and what is left over is systematically flatter than the process it came from.

3. Iterate the trend fit to generalised least squares #

The fix is to feed the fitted covariance back into the trend estimate and go round again. Two or three passes is normally enough.

D = squareform(pdist(coords))

def exp_cov(D, psill, eff_range, nugget):
    """Exponential covariance; eff_range is the 95% (practical) range."""
    C = psill * np.exp(-D / (eff_range / 3.0))
    C[np.diag_indices_from(C)] += nugget
    return C

params, beta = (rng_o, psill_o, nug_o), beta_ols
for it in range(3):
    C = exp_cov(D, params[1], params[0], params[2])
    CiF = np.linalg.solve(C, F)
    beta = np.linalg.solve(F.T @ CiF, CiF.T @ soc)
    Vg = skg.Variogram(coords, soc - F @ beta, model="exponential",
                       n_lags=15, maxlag=6000, normalize=False)
    params = Vg.parameters
    print(f"iter {it}: b0={beta[0]:.4f} b1={beta[1]:.6f}  "
          f"range={params[0]:.0f} psill={params[1]:.4f} nugget={params[2]:.4f}")

iter 0: b0=1.1204 b1=0.008731  range=1352 psill=0.3178 nugget=0.0538
iter 1: b0=1.1590 b1=0.008662  range=1424 psill=0.3341 nugget=0.0523
iter 2: b0=1.1628 b1=0.008655  range=1431 psill=0.3356 nugget=0.0521

The slope moves from 0.009214 to 0.008655 against a truth of 0.0085, and the variogram converges to a range of 1431 m and a total sill of 0.3877 — within 5% and 3% of the truth respectively. The trend coefficient changed by 6%; the variogram range changed by 28%. That asymmetry is the practical message: OLS is a defensible trend estimator and a poor residual generator.

4. Run universal kriging, then reproduce it by hand #

With the converged variogram in place, both routes can be run on identical inputs. gs.krige.ExtDrift is external drift kriging, which is universal kriging with a covariate rather than a coordinate polynomial as the drift term.

g1 = np.linspace(0.0, SIDE, 60)
GX, GY = np.meshgrid(g1, g1)
gx, gy = GX.ravel(), GY.ravel()
elev_g = elevation(gx, gy)

psill, eff_range, nug = 0.3356, 1431.0, 0.0521
model = gs.Exponential(dim=2, var=psill, len_scale=eff_range / 3.0, nugget=nug)

uk = gs.krige.ExtDrift(model, cond_pos=(xs, ys), cond_val=soc,
                       ext_drift=elev_s, exact=True)
uk((gx, gy), mesh_type="unstructured", ext_drift=elev_g)
uk_pred, uk_var = uk.field.copy(), uk.krige_var.copy()

The regression kriging equivalent is the two-line predictor from the algebra above, written out with numpy.

C = exp_cov(D, psill, eff_range, nug)
CiF = np.linalg.solve(C, F)
A = F.T @ CiF                                   # F' C^-1 F
beta_gls = np.linalg.solve(A, CiF.T @ soc)
resid = soc - F @ beta_gls

Dg = np.hypot(gx[:, None] - xs[None, :], gy[:, None] - ys[None, :])
c0 = psill * np.exp(-Dg / (eff_range / 3.0))    # no nugget: s0 is not a sample
f0 = np.column_stack([np.ones(gx.size), elev_g])

rk_pred = f0 @ beta_gls + c0 @ np.linalg.solve(C, resid)

sk_var = (psill + nug) - np.einsum("ij,ij->i", c0, np.linalg.solve(C, c0.T).T)
d0 = f0 - c0 @ CiF
trend_var = np.einsum("ij,ij->i", d0 @ np.linalg.inv(A), d0)

print(f"max |RK(GLS) - UK| prediction = {np.abs(rk_pred - uk_pred).max():.3e}")
print(f"max |RK(GLS) - UK| variance   = "
      f"{np.abs(sk_var + trend_var - uk_var).max():.3e}")

max |RK(GLS) - UK| prediction = 4.441e-13
max |RK(GLS) - UK| variance   = 9.770e-13

That is the equivalence, demonstrated rather than asserted. A cross-check against a third implementation is worth the four extra lines; see universal kriging with external drift in PyKrige for the full treatment of that API.

from pykrige.uk import UniversalKriging

# PyKrige's "sill" is the TOTAL sill, gstools' "var" is the PARTIAL sill.
uk_pk = UniversalKriging(
    xs, ys, soc, variogram_model="exponential",
    variogram_parameters={"sill": psill + nug, "range": eff_range, "nugget": nug},
    drift_terms=["specified"], specified_drift=[elev_s])
z_pk, _ = uk_pk.execute("points", gx, gy, specified_drift_arrays=[elev_g])
print(f"max |UK(pykrige) - UK(gstools)| = {np.abs(z_pk.data - uk_pred).max():.3e}")

max |UK(pykrige) - UK(gstools)| = 2.8e-10

5. Compare the naive regression kriging output against universal kriging #

Now run the version most people actually ship: OLS trend, OLS-residual variogram, simple kriging of the residuals, and the residual kriging variance reported as the prediction variance.

C_o = exp_cov(D, psill_o, rng_o, nug_o)
c0_o = psill_o * np.exp(-Dg / (rng_o / 3.0))
rk_naive = f0 @ beta_ols + c0_o @ np.linalg.solve(C_o, r_ols)
naive_var = (psill_o + nug_o) - np.einsum(
    "ij,ij->i", c0_o, np.linalg.solve(C_o, c0_o.T).T)

print(f"max |RK(OLS) - UK| = {np.abs(rk_naive - uk_pred).max():.3e}")
print(f"rms |RK(OLS) - UK| = {np.sqrt(np.mean((rk_naive - uk_pred)**2)):.3e}")
print(f"mean SK-of-residual variance (same variogram) = {sk_var.mean():.4f}")
print(f"mean UK variance                              = {uk_var.mean():.4f}")
print(f"mean ratio = {(sk_var / uk_var).mean():.3f}   "
      f"worst node = {(sk_var / uk_var).min():.3f}")
print(f"mean naive variance (OLS variogram too) = {naive_var.mean():.4f}")

max |RK(OLS) - UK| = 4.10e-02
rms |RK(OLS) - UK| = 9.20e-03
mean SK-of-residual variance (same variogram) = 0.1663
mean UK variance                              = 0.1885
mean ratio = 0.883   worst node = 0.711
mean naive variance (OLS variogram too) = 0.1449

Interpreting the Output #

The predictions barely move. Across 3600 grid nodes the OLS-trend version differs from universal kriging by at most 0.041 percentage points of organic carbon, against a data range of 1.11 to 5.02. Nobody would notice that on a map, and no cross-validation on 180 points would resolve it. If the only thing you ship is a surface, the choice of route is close to irrelevant.

The variances move a great deal. Holding the variogram fixed so that the trend term is isolated, the simple-kriging residual variance averages 0.1663 against a universal kriging variance of 0.1885 — a mean ratio of 0.883. A 95% prediction interval built on the naive number has a half-width of $1.96\sqrt{0.1663} = 0.80$ where the correct half-width is $1.96\sqrt{0.1885} = 0.85$. That is a 6% shortfall on average, and at the worst node the ratio falls to 0.711, a 16% shortfall. Compound that with the biased OLS-residual variogram and the mean naive variance drops to 0.1449, 23% below the correct value.

What good looks like: a residual variogram whose empirical points sit on a clean rising curve with no upward drift at long lags, GLS trend coefficients that stop moving after two iterations, and a ratio of naive to full variance above about 0.9 across the prediction domain. Warning signs are a residual variogram that keeps climbing at the largest lags, meaning the trend has not removed the drift; trend coefficients that swing by more than their standard errors between iterations, meaning the covariance model is unstable; and a variance ratio below 0.8 anywhere you intend to report an interval.

6. The case where they genuinely diverge #

Change the truth so that organic carbon saturates above 250 m — a common enough soil response — and a drift linear in elevation can no longer represent it. Regression kriging can, because the trend model is arbitrary.

def rk_fold(trend, tr, te, model):
    trend.fit(elev_s[tr, None], soc_nl[tr])
    r = soc_nl[tr] - trend.predict(elev_s[tr, None])
    sk = gs.krige.Simple(model, cond_pos=(xs[tr], ys[tr]), cond_val=r,
                         mean=0.0, exact=True)
    sk((xs[te], ys[te]), mesh_type="unstructured")
    return trend.predict(elev_s[te, None]) + sk.field

def uk_fold(tr, te, model):
    k = gs.krige.ExtDrift(model, cond_pos=(xs[tr], ys[tr]), cond_val=soc_nl[tr],
                          ext_drift=elev_s[tr], exact=True)
    k((xs[te], ys[te]), mesh_type="unstructured", ext_drift=elev_s[te])
    return k.field

soc_nl = 1.20 + 0.0085 * np.minimum(elev_s, 250.0) + e_corr \
         + rng.normal(0.0, np.sqrt(0.05), N)

kf = KFold(n_splits=5, shuffle=True, random_state=0)
for name, fn in [("universal kriging", lambda tr, te: uk_fold(tr, te, model)),
                 ("RK + OLS trend",
                  lambda tr, te: rk_fold(LinearRegression(), tr, te, model)),
                 ("RK + random forest",
                  lambda tr, te: rk_fold(
                      RandomForestRegressor(n_estimators=300, min_samples_leaf=5,
                                            random_state=0), tr, te, model))]:
    err = np.concatenate([fn(tr, te) - soc_nl[te] for tr, te in kf.split(xs)])
    print(f"{name:<20} RMSE = {np.sqrt((err**2).mean()):.3f}")

universal kriging    RMSE = 0.631
RK + OLS trend       RMSE = 0.629
RK + random forest   RMSE = 0.487

Repeat the same loop on the original linear soc and the ordering reverses: universal kriging scores 0.402, the OLS trend 0.404, and the random forest 0.451. The forest pays for its flexibility with variance, and when the truth is a straight line that payment buys nothing. This is the real content of the choice, and it is the same trade-off discussed under combining trend models with kriging residuals.

Critical Best Practices #

Fit the residual variogram to GLS residuals, not OLS residuals #

This is the single correction that changes numbers rather than opinions. Fit once with OLS to get a starting covariance, then loop: build $C$, refit $\hat{\beta}$ by GLS, refit the variogram on the new residuals, repeat until the parameters settle. Three passes cost milliseconds on a few hundred points and recovered a range 28% closer to the truth in the run above. Restricted maximum likelihood is the more principled version of the same idea and is worth the extra machinery when $n$ is small.

Never report the simple-kriging residual variance as the prediction variance #

If you take the regression kriging route and then quote the residual kriging variance, your intervals are too narrow everywhere, and worst exactly where users are most likely to be misled — at the edges of the survey and where the covariate is extrapolated. Either add the $d_0^{\mathsf{T}} (F^{\mathsf{T}} C^{-1} F)^{-1} d_0$ term as shown, or, if the trend model is non-linear and no closed form exists, bootstrap the whole two-stage fit and take empirical quantiles.

The words “range”, “sill” and “len_scale” are not portable between libraries #

gstools parameterises the exponential model with len_scale, the e-folding scale, so the practical range is 3 * len_scale, and its var is the partial sill above the nugget. pykrige and scikit-gstat both take the effective range directly, but pykrige’s sill key means the total sill and it derives the partial sill as sill - nugget. Passing a total sill where a partial sill is expected inflates the variance by exactly the nugget and is invisible in the predictions, which is why it survives code review.

Check the covariate is available and in range everywhere you predict #

External drift kriging needs the drift value at every prediction node, not just at the samples, and the unbiasedness constraints extrapolate the trend without complaint. If elevation on the grid runs from 80 m to 338 m but the samples only cover 113–330 m, the trend is being used outside its support at both ends. The trend-uncertainty term is what tells you this is happening — it is the term that grew to 29% of the total at the extrapolating node in the chart above — which is another argument for computing it even when you take the regression kriging route.

Cross-validate the trend model separately from the kriging #

A flexible trend model can fit the residual structure rather than the trend, leaving a residual variogram that is pure nugget and a kriging step that does nothing. Fit the trend inside the cross-validation folds, as the code above does, never once on the full dataset. If the residual variogram of an out-of-fold trend is flat while the in-sample one is not, the trend has eaten the spatial structure and the split between deterministic and stochastic components needs revisiting.

Troubleshooting #

Next Steps #

If the trend needs a model class kriging cannot express, work through combining trend models with kriging residuals for the full two-stage pipeline with honest uncertainty; if the drift is linear, go straight to universal kriging with external drift in PyKrige.

Frequently Asked Questions #

Are regression kriging and universal kriging the same thing?

They are the same estimator whenever the trend is linear in the drift terms, the trend coefficients come from generalised least squares, and both use the same covariance model. Under those conditions the predictions agree to machine precision, because both are the best linear unbiased predictor for the same model. Regression kriging as usually implemented differs only because the trend is fitted by ordinary least squares instead, and because the reported variance omits the term covering uncertainty in the trend coefficients.

Why is a variogram fitted to OLS residuals biased?

Ordinary least squares minimises the residual sum of squares while ignoring the correlation between residuals, so it absorbs part of the long-range spatially correlated variation into the fitted trend. What is left over is systematically too flat: the empirical variogram of OLS residuals rises to a lower sill and reaches it sooner. In the worked example the sill came out eighteen per cent low and the range twenty-five per cent short. Iterating the trend fit to generalised least squares recovers both to within five per cent.

Does regression kriging really understate the prediction variance?

The naive implementation does, because it reports the kriging variance of the residual field alone and treats the fitted trend as if it were known exactly. Universal kriging adds a second term that grows with distance from the data and with extrapolation of the covariate. Across the example grid the naive variance averaged eighty-eight per cent of the universal kriging variance, falling to seventy-one per cent where the elevation covariate ran beyond its sampled range. Prediction intervals shrink by roughly six per cent on average and sixteen per cent at worst.

When should I choose regression kriging over universal kriging?

Choose regression kriging when the trend cannot be written as a linear combination of drift terms — a random forest, a gradient boosting model, a generalised additive model with smooths — or when you need to inspect and model the residual variogram by eye. Choose universal kriging when the drift is linear in its coefficients and you need a defensible prediction variance, because one system delivers both without bookkeeping. If the trend is linear the two give the same answer, so the choice is operational.

Related

• Combining Trend Models with Kriging Residuals — the two-stage pipeline in full, including non-linear trend models
• Universal Kriging with External Drift in PyKrige — the single-system route with a covariate drift
• Stationarity & Trend Analysis — deciding what counts as trend before either method is fitted

← Back to Regression Kriging