# Poor performance from gstat::krige with a noise predictor

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  ## 1 Answer 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")  # 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")


bubble(pred, "resid")


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")


bubble(pred, "resid")


sqrt(mean((pred\$resid)^2))

## [4] 2.171701