# How can I create a map of variance for IDW predictions?

Unlike with kriging, predicting spatial variation using the IDW function from gstat returns only predicted values, not estimates of variation.

library(sp)
data("meuse")
coordinates(meuse) = ~x+y
data("meuse.grid")
gridded(meuse.grid) = ~x+y
m.idw <- gstat::idw(zinc~1, meuse, newdata=meuse.grid, idp=2)


There is a column for variation, but it is empty. Is this a fundamental mathematical limitation of the IDW approach, or is it just not implemented in gstat, or is there something else wrong in my approach?

• Yes, yes, and impossible to tell. IDW is non-statistical: it's just an algorithm; the weights it uses do not depend on the data values at all. What you can do is use a variogram along with the IDW weights to estimate a prediction variance--which can be a useful way to see what using the simpler, faster IDW algorithm might be costing you compared to kriging. – whuber Jul 30 '19 at 21:38
• Thanks @whuber, this is exactly what I needed to know. If you make it an answer I'll accept. – Nat Jul 30 '19 at 22:59

I will supply two answers: first, this is a limitation of IDW; but second, there are ways to overcome it by making assumptions.

IDW (inverse distance weighting) is a procedure--best used in an exploratory mode--for predicting unobserved values of a function defined on a metric space such as the Euclidean plane or sphere. Specifically, suppose you have observed $$y_i = f(x_i)$$ at a set of distinct points $$x_i\in \mathcal{X}$$ where $$\mathcal X$$ has a metric $$\delta$$ and you wish to predict the values of $$f$$ at other points $$x^\prime_j\in \mathcal X$$ (which often are chosen to support the creation of graphs, contour maps, images, etc. of $$f$$).

Carrying out IDW requires two choices of the analyst (both of which are often made arbitrarily by accepting software defaults): (1) selecting a power $$p\ge 1$$ and (2) defining a neighborhood search procedure. The latter is a way of selecting observation points near any given point $$x.$$ For our purposes the details don't matter, so let's just write $$\mathcal{N}(x) \subset \{x_1,x_2,\ldots, x_n\}$$ for the neighborhood of $$x\in\mathcal{M}.$$

Given these choices, IDW is just a recipe: from the distances $$\delta(x,x_i)$$ between $$x$$ and its neighboring observations $$x_i\in\mathcal{N}(x),$$ construct relative weights

$$w_i(x) = \delta(x,x_i)^{-p} \gt 0,$$

normalize them to sum to unity by setting

$$\omega_i(x) = \frac{w_i(x)}{\sum_{x_i\in\mathcal{N}(x)} w_i(x)},$$

and predict the function's value as a weighted combination of the neighboring observations,

$$\hat y(x) =\sum_{x_i\in\mathcal{N}(x)} \omega_i(x)\, y_i.$$

Because there is no probability model involved, there is no general way to estimate a variance or standard error of estimate.

One could use various cross-validation tactics to estimate prediction errors. These all implicitly assume $$f$$ has various restrictive properties concerning how rapidly it can change from point to point. Rather than dodging considerations of what is actually being assumed about $$f,$$ I have found it useful at times to do a quick-and-dirty geostatistical analysis of the data, which (although it also makes assumptions) can reveal something about the properties of $$f.$$

This means supposing $$f$$ is a realization of some kind of stochastic process $$F$$ ("random field") on $$\mathcal{M},$$ so that $$(y_1,y_2,\ldots,y_n)$$ is a single observation of the $$n$$-variate random variable $$(F(x_1), F(x_2), \ldots, F(x_n));$$ and estimating its covariance function

$$C(x, x^\prime) = \operatorname{Cov}(F(x), F(x^\prime)).$$

This is usually accomplished by adopting some assumptions about stationarity and ergodicity so that the data can be used to estimate $$C$$ (in the usual way, likely through some form of variogram estimation). At other times, you might use estimates of $$C$$ based on related data.

With such a covariance function in hand it's a relatively simple matter to compute the variance of any predicted value as $$\operatorname{Var}(\hat{y}(x)) = \operatorname{Var}\left(\sum_{x_i\in\mathcal{N}(x)} \omega_i(x)\, y_i\right) = \sum_{i,j\in\mathcal{N}(x)} \omega_i(x)\omega_j(x)\,C(x_i,x_j).$$

Its square root serves as an estimated standard error of prediction.

Note that what makes the foregoing formula possible is that the weights $$\omega_i(x)$$ depend only on the locations of the observations, which are considered determinate (non-random), and not on the observations themselves.

Among other applications, this approach affords a principled way to adjust the otherwise arbitrary choices of power $$p$$ and neighborhood search $$\mathcal{N}$$ made at the outset: if one set of choices consistently results in smaller standard errors (in regions of $$\mathcal M$$ of particular interest) than another set, then the former is to be preferred. Note that making such a comparison does not require excellent estimation of $$C$$; in particular, the hardest part of estimating $$C$$ usually is finding an appropriate scale factor--but that's irrelevant if one is only comparing two predictions based on the same $$C.$$ Usually, determining the right range of the covariance function is what matters most. That's why this approach usually is not onerous and can produce useful information.

As others already pointed out, IDW is not based on a theory that could provide variance maps. On the other hand, it is easy to arbitrarily choose a metric, and use it to quantify how far from the observations the estimates are. This is what I have implemented in GeoStats.jl:

using GeoStats
using InverseDistanceWeighting
using Plots

sdata   = PointSetData(Dict(:z => [1.,0.,1.]), [25. 50. 75.; 25. 75. 50.])
sdomain = RegularGrid{Floa64}(100, 100)

problem = EstimationProblem(sdata, sdomain, :z)

solver = InvDistWeight()

solution = solve(problem, solver)

contourf(solution)


And this is how it looks like if you choose a Chebyshev metric:

solver = InvDistWeight(:z => (distance=Chebyshev(),))

solution = solve(problem, solver)

contourf(solution)


• Because this is interesting, could you describe or explain what the code does and what the "variable variance" maps actually represent? – whuber Aug 5 '19 at 13:51
• Hi @whuber, you can find the full notebook with the example here github.com/juliohm/InverseDistanceWeighting.jl/blob/master/docs/…. The name variable is copied from the notebook instead, sorry for the confusion. The only reason I called the title of the plot variable variance is because I couldn't find a better, more general, name for it. In Kriging we get variance, but for other geostatistical methods implemented in GeoStats.jl we can still get some "uncertainty" map that reflects how far we are from the observations. – juliohm Aug 6 '19 at 11:49
• Some demonstration that these maps actually do reflect some quantitative notion of uncertainty would be welcome. – whuber Aug 6 '19 at 12:02
• In my comment I wrote "uncertainty" in between quotes for a reason. There is no formal definition of uncertainty in this plot. The maps are just to indicate how far from the observations the estimates are according to a metric. In other solvers like locally weighted regression for example, we get a more formal definition, which is derived analytically from the normal equations. Since inverse distance weighting is a heuristic, we can also use a heuristic for the "uncertainty" map. If you have another heuristic in mind that would be more useful, please share. – juliohm Aug 7 '19 at 21:32
• I already did in my reply to this question. – whuber Aug 7 '19 at 23:05