1
$\begingroup$

I'm new to kriging, and I'm considering replacing a use of inverse-distance weighting (IDW) in a spatial modeling project (implemented with gstat::idw in R) with a kriging-based approach (implemented with gstat::krige). On my first attempt, I got poor predictive accuracy. I've been playing around with simulations to figure out what's going on, and it looks like when krige gets one noise predictor, compared to no predictors, it tanks. Is this a bug in gstat or some way in which universal kriging is vastly more prone to overfitting than, say, linear regression? What's going on?

Here's an example:

library(gstat)

set.seed(1000)
n = 1000
d = data.frame(
    train = c(rep(T, .7*n), rep(F, .3*n)),
    x = runif(n, 0, 10),
    y = runif(n, 0, 10),
    noise = rnorm(n))
d = transform(d, outcome = (x - 5)^2 + (y - 5)^2 + rnorm(n, sd = 2))

obs = d[!d$train, "outcome"]

tryit = function(my.formula)
   {preds = gstat::krige(
        formula = my.formula,
        locations = ~ x + y,
        data = d[d$train,],
        newdata = d[!d$train,],
        debug.level = 0)[, "var1.pred"]
    message("RMSE: ", sqrt(mean((obs - preds)^2)))}

message("SD: ", sd(obs))
tryit(outcome ~ 1)
tryit(outcome ~ 1 + noise)

The output is:

SD: 10.8120249379363
RMSE: 3.68179409633032
RMSE: 10.7943966111937
$\endgroup$
1
$\begingroup$

I think the problem was lack of a variogram model. I'm not actually sure how krige handles this situation, but perhaps it has a default variogram model fitting procedure which is inadequate.

In any case, a variogram model can be automatically fitted using package automap, which is the easiest option.

Please see below two calculation of RMSE - 1. Ordinary Kriging (outcome ~ 1) 2. Universal Kriging (outcome ~ noise)

As could be expected, the RMSE is almost identical in both cases (RMSE=2.17) since the noise has no contribution in explaining the trend.

library(sp)
library(gstat)
library(automap)

set.seed(1000)
n = 1000
d = data.frame(
    train = c(rep(T, .7*n), rep(F, .3*n)),
    x = runif(n, 0, 10),
    y = runif(n, 0, 10),
    noise = rnorm(n))
d = transform(d, outcome = (x - 5)^2 + (y - 5)^2 + rnorm(n, sd = 2))

# To points
pnt = d
coordinates(pnt) = ~ x + y

# Observed
bubble(pnt, "outcome")

enter image description here

# Ordinary Kriging
f = outcome ~ 1
v = autofitVariogram(formula = f, input_data = pnt)
g = gstat(formula = f, data = pnt[pnt$train, ], model = v$var_model)
pred = predict(object = g, newdata = pnt[!pnt$train, ])
## [using ordinary kriging]
pred$resid = pnt$outcome[!pnt$train] - pred$var1.pred
bubble(pred, "var1.pred")

enter image description here

bubble(pred, "resid")

enter image description here

sqrt(mean((pred$resid)^2))
## [4] 2.173127
# Universal Kriging
f = outcome ~ noise
v = autofitVariogram(formula = f, input_data = pnt)
g = gstat(formula = f, data = pnt[pnt$train, ], model = v$var_model)
pred = predict(object = g, newdata = pnt[!pnt$train, ])
## [using universal kriging]
pred$resid = pnt$outcome[!pnt$train] - pred$var1.pred
bubble(pred, "var1.pred")

enter image description here

bubble(pred, "resid")

enter image description here

sqrt(mean((pred$resid)^2))
## [4] 2.171701
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.