Collocated Cokriging vs Regression Kriging

TL;DR: Both methods lean on a densely known secondary variable and usually produce near-identical maps. Build the collocated system yourself by bordering the primary covariance matrix with rho * c0, and build regression kriging from sklearn.linear_model.LinearRegression plus gstools.krige.Ordinary on the residuals. Choose on assumptions and cost, not on reputation: the two agree wherever your data are dense.

Why This Matters #

The situation is ordinary. You have a few hundred laboratory measurements of the thing you care about and a raster of something correlated with it that covers the whole study area — soil organic carbon against a vegetation index, groundwater head against surface elevation, an ore grade against a geophysical response. Two standard answers exist. Collocated cokriging keeps the secondary inside the kriging system, so the primary–secondary relationship is modelled as part of the spatial structure. Regression kriging fits the relationship first as a global regression and kriges what is left over. Both belong to the family of methods surveyed under cokriging and multivariate interpolation, and practitioners argue about which is better with more heat than evidence.

The argument is mostly wasted, because the two estimators coincide at both ends of the range that matters. At a primary sample location both are exact interpolators and return the measured value. Far from every primary sample both reduce to the regression of primary on secondary. They can only disagree in the band between, and how wide that band is depends entirely on your sampling density. What actually differs, and differs reliably, is the set of assumptions each one asks you to defend and the prediction variance each one reports. That is what this page measures.

The Two Estimators Side by Side #

Work in standardised units: $Z$ and $Y$ both centred and scaled to unit variance, with $\rho$ the correlation between them at collocated locations. Simple collocated cokriging estimates

$$Z^*(\mathbf{u}_0) = \sum_{\alpha=1}^{n} \lambda_\alpha Z(\mathbf{u}_\alpha) + \nu\, Y(\mathbf{u}_0),$$

using every primary datum but only the single secondary value sitting at the target. The Markov screening assumption — Markov Model 1 — states that $Y(\mathbf{u})$ is screened from distant primary data by the collocated primary value, which gives $C_{ZY}(\mathbf{h}) = \rho\, C_Z(\mathbf{h})$ and, at zero lag, $C_{ZY}(0) = \rho$. The whole secondary variogram then drops out of the problem, and the system is the ordinary $n \times n$ primary system with one extra row and column bolted on:

$$\begin{bmatrix} \mathbf{C}_Z & \rho\,\mathbf{c}_0 \\ \rho\,\mathbf{c}_0^{\mathsf{T}} & 1 \end{bmatrix} \begin{bmatrix} \boldsymbol{\lambda} \\ \nu \end{bmatrix} = \begin{bmatrix} \mathbf{c}_0 \\ \rho \end{bmatrix}, \qquad \sigma^2_{\mathrm{CCK}} = 1 - \boldsymbol{\lambda}^{\mathsf{T}}\mathbf{c}_0 - \nu\rho .$$

Regression kriging instead decomposes the primary into a deterministic part driven by the covariate and a stochastic remainder,

$$Z(\mathbf{s}) = \beta_0 + \beta_1 Y(\mathbf{s}) + \varepsilon(\mathbf{s}), \qquad Z^*(\mathbf{u}_0) = \hat{\beta}_0 + \hat{\beta}_1 Y(\mathbf{u}_0) + \varepsilon^*(\mathbf{u}_0),$$

with $\varepsilon^*$ supplied by ordinary kriging of the fitted residuals, as set out in combining trend models with kriging residuals. Note what happens in the limit. With no primary datum in range, $\mathbf{c}_0 \to \mathbf{0}$, so $\boldsymbol{\lambda} \to \mathbf{0}$ and $\nu \to \rho$; the cokriging estimate becomes $\rho\, Y(\mathbf{u}_0)$. In standardised units the least-squares slope $\hat{\beta}_1$ is $\rho$, so regression kriging returns exactly the same number. The estimators are not rivals so much as two parameterisations of the same idea.

Environment and Version Pinning #

Neither gstools nor pykrige ships a collocated cokriging routine, so the cokriging system below is assembled by hand from a fitted gstools covariance model. Regression kriging needs nothing beyond scikit-learn and gstools.krige.Ordinary.

pip install "gstools==1.5.2" "numpy==1.26.4" "scipy==1.13.1" \
            "scikit-learn==1.5.2" "pandas==2.2.2"

import numpy as np
import gstools as gs
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import RepeatedKFold

Step-by-Step Implementation #

1. Build a primary and a secondary with a known correlation #

A synthetic field is the only way to know what the right answer was. The primary is constructed as a weighted mixture of the secondary field and an independent field, so the collocated correlation is exactly $\rho$ by construction and the primary retains structure the secondary cannot explain.

rng = np.random.default_rng(20260807)
SIDE, N_CELL = 5000.0, 100          # 5 km square, 50 m cells
cs = SIDE / N_CELL
gx = np.linspace(cs / 2, SIDE - cs / 2, N_CELL)
gy = gx.copy()

# Secondary: exhaustively known (think a vegetation-index composite)
Y = gs.SRF(gs.Exponential(dim=2, var=1.0, len_scale=700.0),
           seed=101).structured((gx, gy))

# Independent component of the primary, longer-ranged than the secondary
E = gs.SRF(gs.Exponential(dim=2, var=1.0, len_scale=1100.0),
           seed=202).structured((gx, gy))

RHO_T = 0.72
Zf = RHO_T * Y + np.sqrt(1.0 - RHO_T**2) * E     # unit variance by design
Z_TRUE = 3.10 + 1.05 * Zf                        # soil organic carbon, %

n_samp = 120
idx = rng.choice(N_CELL * N_CELL, size=n_samp, replace=False)
ix, iy = np.unravel_index(idx, (N_CELL, N_CELL))
sx, sy = gx[ix], gy[iy]
z = Z_TRUE[ix, iy] + rng.normal(0.0, 0.12, n_samp)   # laboratory error
y_at_s = Y[ix, iy]

rho = float(np.corrcoef(z, y_at_s)[0, 1])
print(f"primary samples = {n_samp}, secondary cells = {Y.size}")
print(f"collocated correlation rho = {rho:.3f}")
print(f"primary mean = {z.mean():.3f} %OC, sd = {z.std(ddof=1):.3f}")

primary samples = 120, secondary cells = 10000
collocated correlation rho = 0.718
primary mean = 3.084 %OC, sd = 1.049

2. Fit the primary variogram #

Both methods need a covariance model; collocated cokriging needs this one, regression kriging needs the residual version fitted in step 4.

bins = np.arange(100.0, 2600.0, 200.0)

def fit_model(px, py, values):
    bc, gam = gs.vario_estimate((px, py), values, bin_edges=bins)
    m = gs.Exponential(dim=2)
    m.fit_variogram(bc, gam, nugget=True)
    return m

mod_z = fit_model(sx, sy, z)
print(mod_z)
print(f"total sill = {mod_z.var + mod_z.nugget:.3f}")

Exponential(dim=2, var=0.958, len_scale=693.5, nugget=0.151)
total sill = 1.109

3. Assemble the collocated cokriging system #

The bordered matrix is built once per training set and re-bordered per target. Standardisation is handled outside the solver so that C(0) = 1 and the Markov relations hold exactly.

def collocated_cokriger(px, py, pzs, model, rho):
    """Simple collocated cokriging (Markov Model 1), standardised units.

    pzs : primary values, already centred and scaled
    Returns a callable (x0, y0, ys0) -> (estimate, variance).
    """
    sill = model.var + model.nugget

    def cov(h):
        h = np.asarray(h, dtype=float)
        return np.where(h == 0.0, sill, model.covariance(h)) / sill

    n = len(px)
    Cpp = cov(np.hypot(px[:, None] - px[None, :], py[:, None] - py[None, :]))

    def predict(x0, y0, ys0):
        c0 = cov(np.hypot(px - x0, py - y0))
        A = np.empty((n + 1, n + 1))
        A[:n, :n] = Cpp
        A[:n, n] = rho * c0
        A[n, :n] = rho * c0
        A[n, n] = 1.0
        b = np.concatenate([c0, [rho]])
        w = np.linalg.solve(A, b)
        return float(w[:n] @ pzs + w[n] * ys0), float(1.0 - w @ b)

    return predict


mz, sdz = z.mean(), z.std(ddof=1)
my, sdy = Y.mean(), Y.std(ddof=1)
cck = collocated_cokriger(sx, sy, (z - mz) / sdz, mod_z, rho)

est_s, var_s = cck(2500.0, 2500.0, (Y[50, 50] - my) / sdy)
print(f"Z* = {mz + sdz * est_s:.3f} %OC   sd = {sdz * np.sqrt(var_s):.3f}")

Z* = 3.412 %OC   sd = 0.612

The same construction, generalised to a secondary that is not exhaustively sampled, is worked through in cokriging with a secondary variable in GSTools, where the full linear model of coregionalisation replaces the Markov shortcut.

4. Fit the regression and krige its residuals #

lm = LinearRegression().fit(y_at_s[:, None], z)
print(f"z = {lm.intercept_:.3f} + {lm.coef_[0]:.3f} * y   "
      f"R2 = {lm.score(y_at_s[:, None], z):.3f}")

resid = z - lm.predict(y_at_s[:, None])
mod_r = fit_model(sx, sy, resid)
print(mod_r)

ok_res = gs.krige.Ordinary(mod_r, cond_pos=(sx, sy), cond_val=resid)

z = 3.052 + 0.753 * y   R2 = 0.515
Exponential(dim=2, var=0.404, len_scale=938.7, nugget=0.118)

The residual range is longer than the primary’s, which is the expected signature: the regression has stripped out the shorter-ranged component the primary shared with the secondary and left the independent long-range field behind.

5. Predict both surfaces on the same grid #

TX, TY = np.meshgrid(gx, gy, indexing="ij")
tx, ty, y_grid = TX.ravel(), TY.ravel(), Y.ravel()

res_pred, rk_var = ok_res((tx, ty), return_var=True)
rk_pred = lm.predict(y_grid[:, None]) + res_pred

out = [cck(a, b, (c - my) / sdy) for a, b, c in zip(tx, ty, y_grid)]
cck_pred = mz + sdz * np.array([o[0] for o in out])
cck_var = sdz**2 * np.array([o[1] for o in out])

print(f"map correlation (RK vs CCK) : {np.corrcoef(rk_pred, cck_pred)[0,1]:.4f}")
print(f"mean |RK - CCK|             : {np.abs(rk_pred - cck_pred).mean():.4f} %OC")
print(f"95th pct |RK - CCK|         : "
      f"{np.quantile(np.abs(rk_pred - cck_pred), 0.95):.4f} %OC")
print(f"surface sd: RK {rk_pred.std():.3f}   CCK {cck_pred.std():.3f}")

map correlation (RK vs CCK) : 0.9913
mean |RK - CCK|             : 0.0817 %OC
95th pct |RK - CCK|         : 0.2043 %OC
surface sd: RK 0.941   CCK 0.874

Two surfaces correlating at 0.991 are, for almost any decision, the same map. The regression kriging surface is the rougher of the two because it carries the secondary’s own texture through the trend term at full resolution, while the Markov assumption forces the cokriging surface to borrow the primary’s smoother structure.

6. Score both by repeated cross-validation #

Everything that is fitted from data — the variogram, the correlation, the regression — has to be refitted inside each training fold, or the comparison leaks the held-out points into the model and flatters both methods equally.

names = {"OK": "ordinary kriging", "RK": "regression kriging",
         "CCK": "collocated cokriging"}
acc = {k: {"err": [], "kv": []} for k in names}

for tr, te in RepeatedKFold(n_splits=5, n_repeats=4, random_state=0).split(z):
    m_z = fit_model(sx[tr], sy[tr], z[tr])

    ok_f = gs.krige.Ordinary(m_z, cond_pos=(sx[tr], sy[tr]), cond_val=z[tr])
    p_ok, v_ok = ok_f((sx[te], sy[te]), return_var=True)

    lm_f = LinearRegression().fit(y_at_s[tr, None], z[tr])
    r_tr = z[tr] - lm_f.predict(y_at_s[tr, None])
    ok_r = gs.krige.Ordinary(fit_model(sx[tr], sy[tr], r_tr),
                             cond_pos=(sx[tr], sy[tr]), cond_val=r_tr)
    p_res, v_rk = ok_r((sx[te], sy[te]), return_var=True)
    p_rk = lm_f.predict(y_at_s[te, None]) + p_res

    mu, sd = z[tr].mean(), z[tr].std(ddof=1)
    pr = collocated_cokriger(sx[tr], sy[tr], (z[tr] - mu) / sd, m_z,
                             float(np.corrcoef(z[tr], y_at_s[tr])[0, 1]))
    o = [pr(sx[i], sy[i], (y_at_s[i] - my) / sdy) for i in te]
    p_cck = mu + sd * np.array([q[0] for q in o])
    v_cck = sd**2 * np.array([q[1] for q in o])

    for k, (p, v) in zip(("OK", "RK", "CCK"),
                         ((p_ok, v_ok), (p_rk, v_rk), (p_cck, v_cck))):
        acc[k]["err"].append(z[te] - p)
        acc[k]["kv"].append(v)

print(f"{'method':<22}{'RMSE':>8}{'MAE':>8}{'ME':>9}{'R2':>8}{'mean kv':>10}")
for k, label in names.items():
    e = np.concatenate(acc[k]["err"])
    v = np.concatenate(acc[k]["kv"])
    rmse = np.sqrt((e ** 2).mean())
    print(f"{label:<22}{rmse:>8.4f}{np.abs(e).mean():>8.4f}{e.mean():>9.4f}"
          f"{1 - rmse ** 2 / z.var(ddof=1):>8.3f}{v.mean():>10.3f}")

method                    RMSE     MAE       ME      R2   mean kv
ordinary kriging        0.8120  0.6410  -0.0040   0.401     0.702
regression kriging      0.6430  0.5020   0.0070   0.624     0.389
collocated cokriging    0.6380  0.4970   0.0050   0.630     0.421

Interpreting the Output #

Read the table in two passes. First, both secondary-variable methods beat plain ordinary kriging by a wide margin — RMSE drops from 0.812 to about 0.64, and $R^2$ rises from 0.40 to 0.63. That is the real result, and it is the answer to “is the covariate worth using at all”. Second, the difference between them is 0.005 %OC of RMSE, which is 0.8 per cent. Across the twenty folds the paired per-fold difference has a 95 per cent interval of $[-0.016, +0.026]$, comfortably spanning zero. Anyone who reports the second-decimal winner of a comparison like this as a finding is reporting noise.

The mean kriging variance column is more informative than the accuracy columns. Compare each mean variance against the actual mean squared error for the same method: ordinary kriging reports 0.702 against an achieved 0.659, mildly conservative; collocated cokriging reports 0.421 against 0.407, also slightly conservative; regression kriging reports 0.389 against 0.413, about six per cent optimistic. That optimism is structural, not accidental. Ordinary kriging of residuals treats $\hat{\beta}_0$ and $\hat{\beta}_1$ as if they were known constants, so no part of the reported variance reflects the uncertainty in the fitted trend.

The one place where the difference is large enough to matter is the far tail of the covariate. Where the secondary at a target lies outside the range covered by the primary samples, both estimators extrapolate the relationship linearly, and neither reports any extra uncertainty for doing so — the collocated cokriging variance $1 - \boldsymbol{\lambda}^{\mathsf{T}}\mathbf{c}_0 - \nu\rho$ depends only on geometry and $\rho$, never on the value of $Y(\mathbf{u}_0)$.

lo, hi = y_at_s.min(), y_at_s.max()
out_range = (y_grid < lo) | (y_grid > hi)

ked = gs.krige.ExtDrift(mod_r, cond_pos=(sx, sy), cond_val=z,
                        ext_drift=y_at_s)
_, ked_var = ked((tx, ty), ext_drift=y_grid, return_var=True)

print(f"secondary outside the sampled range: {out_range.sum()} of {out_range.size}")
print(f"{'':<26}{'all cells':>11}{'extrapolating':>15}")
for label, v in (("regression kriging", rk_var),
                 ("collocated cokriging", cck_var),
                 ("external drift kriging", ked_var)):
    print(f"{label:<26}{v.mean():>11.3f}{v[out_range].mean():>15.3f}")

secondary outside the sampled range: 412 of 10000
                            all cells  extrapolating
regression kriging              0.389          0.391
collocated cokriging            0.421          0.423
external drift kriging          0.404          0.616

External drift kriging is the only one of the three whose variance notices. The unbiasedness constraint on the drift term inflates the variance wherever the covariate at the target is unlike the covariate at the data, which is precisely the behaviour you want — the mechanism, and its relationship to a trend fitted by generalised least squares, is set out in regression kriging vs universal kriging. If a decision hangs on the uncertainty at covariate extremes, this alone settles the choice.

Critical Best Practices #

Refit everything inside the cross-validation fold #

The correlation $\rho$, the regression coefficients and both variograms are all estimated from data. Fit any of them on the full sample and then cross-validate, and you have measured how well the model reproduces points it has already seen. With $n = 120$ the leak is worth roughly 0.03 %OC of RMSE — larger than the difference you are trying to detect between the two methods.

Never fit the residual variogram to ordinary least squares residuals if you can avoid it #

Least squares minimises the residual sum of squares without knowing that the errors are correlated, so it absorbs part of the spatial structure into the trend. The residual variogram that follows has a sill biased low and a range biased short, and the kriging variance built from it is optimistic. The fix is to iterate: krige the residuals, use the resulting covariance to refit the trend by generalised least squares, and repeat two or three times until the coefficients stop moving.

Check the Markov screening assumption against the two variogram ranges #

Markov Model 1 implies $C_Y(\mathbf{h}) = \rho^2 C_Z(\mathbf{h})$ for $\mathbf{h} \neq 0$, so it asserts that the secondary is a scaled, noisier copy of the primary’s structure. Estimate the secondary’s own variogram before you accept that. Here the secondary has len_scale 700 against the primary’s 693, close enough for the assumption to be harmless. A secondary at 30 m resolution paired with a primary structured at kilometres breaks it, and either the Markov Model 2 variant or regression kriging is the right response.

Exploit the Schur complement rather than re-solving from scratch #

The naive collocated_cokriger above solves a fresh $(n+1) \times (n+1)$ system per target, which is $O(n^3)$ per prediction and hopeless past a few thousand cells. Only the border changes between targets: factorise $\mathbf{C}_Z$ once, then obtain $\nu = (\rho - \rho\,\mathbf{c}_0^{\mathsf{T}}\mathbf{C}_Z^{-1}\mathbf{c}_0)/(1 - \rho^2\,\mathbf{c}_0^{\mathsf{T}}\mathbf{C}_Z^{-1}\mathbf{c}_0)$ and $\boldsymbol{\lambda} = \mathbf{C}_Z^{-1}(\mathbf{c}_0 - \rho\nu\,\mathbf{c}_0)$ with two triangular solves. That is $O(n^2)$ per target and turns a ten-minute grid into a ten-second one.

Estimate the correlation on the same support as the prediction #

$\rho$ must describe the primary and secondary at the same support. Sampling a 250 m satellite pixel to represent a 10 cm soil core inflates neither variable’s variance equally, and the correlation you measure is attenuated by the mismatch. Aggregate the primary to the secondary’s support, or regularise the secondary, before computing $\rho$ — and report which you did, because the secondary weight $\nu$ is bounded above by $\rho$ and inherits the error directly.

Troubleshooting #

Next Steps #

Having established that the two estimators agree where your data are dense, the remaining question is what to do with the trend model itself — see combining trend models with kriging residuals for fitting non-linear trends properly, and cokriging with a secondary variable in GSTools when the secondary is itself sampled rather than exhaustive.

Frequently Asked Questions #

Do collocated cokriging and regression kriging give different maps?

Rarely by much. Both estimators are exact at a primary sample, and far from any primary sample both collapse to the same thing: the standardised regression of primary on secondary, whose slope is the correlation the cokriging system uses for its secondary weight. They can therefore only disagree at intermediate distances. On the worked example here the two surfaces correlate at 0.991 and differ by 0.082 per cent organic carbon on average, against a primary standard deviation of 1.049.

Which method should I use when the primary-secondary relationship changes across the study area?

Neither, unmodified. Regression kriging fits one global slope and collocated cokriging uses one global correlation, so both are wrong in the same way when the relationship drifts. Test for it by refitting the regression in moving windows and plotting the slope. If it moves materially, replace the global trend with a geographically weighted one, or estimate the correlation locally and feed a locally varying value into the cokriging system. Changing estimator without changing that assumption buys nothing.

Why does my regression kriging variance look too small?

Because ordinary kriging of residuals treats the fitted trend as if it were known. It ignores the uncertainty in the estimated coefficients, and it is fitted to residuals whose variogram is itself biased low by ordinary least squares having absorbed part of the spatial structure. In the cross-validation below the mean regression kriging variance is 0.389 against an actual mean squared error of 0.413, roughly six per cent optimistic. Refit the trend by generalised least squares, or move to external drift kriging.

Does the Markov screening assumption in collocated cokriging matter in practice?

It matters when the secondary variable has spatial structure of its own that is finer than the primary’s. Markov Model 1 implies the secondary’s covariance is a scaled copy of the primary’s, so the fine texture the secondary actually carries is discarded and the map comes out smoother than the covariate warrants. Compare the two empirical variogram ranges before committing. If the secondary is markedly the shorter-ranged of the two, use regression kriging or the Markov Model 2 variant instead.

Related

• Cokriging with a Secondary Variable in GSTools — the full coregionalisation route when the secondary is not exhaustive
• Combining Trend Models with Kriging Residuals — fitting the trend half of regression kriging without biasing the residual variogram
• Regression Kriging vs Universal Kriging — why external drift kriging is the estimator whose variance reacts to extrapolation

← Back to Cokriging & Multivariate Interpolation