[R] Point estimate from loess contour plot

[Gavin Simpson]

On Thu, 2007-02-08 at 19:09 +0000, Laura Quinn wrote:
> Hi,
> 
> I was wondering if anyone knows of a way by which one can estimate values
> from a contour plot created by using the loess function? I am hoping to
> use the loess contour plot as a means of interpolation to identify
> the loess created values at points at pre-defined (x,y) locations.
> 
> Could anyone point me in the right direction please?
> 
> Thanks.
> 
> Laura Quinn

Hi Laura,

Using an example from MASS (the book by Venables and Ripley, page 423 in
4th Ed (2002)) and the topo data set:

require(MASS)
## loess surface
topo.lo <- loess(z ~ x * y, topo, degree = 1, span = 0.25,
                 normalize = FALSE)
topo.mar <- list(x = seq(0, 6.5, 0.1), y = seq(0, 6.5, 0.1))
new.dat <- expand.grid(topo.mar)
topo.pred <- predict(topo.lo, new.dat)
## draw the contour map based on loess predictions
eqscplot(topo.mar, type = "n")
contour(topo.mar$x, topo.mar$y, topo.pred,
        levels = seq(700, 1000, 25), add = TRUE)
## original points
points(topo, col = "red")

So now we have a loess surface defined by the model (topo.lo) and we can
use the predict.loess method to get point predictions based on the
model. This is what was used to produce the points draw the contour
surface, but on a regular grid. For some new point, not on this regular
grid we can use the same approach:

> predict(topo.lo, data.frame(x = 4.8, y = 3.1))
[1] 824.0046

which yields a height of 824 and a bit feet for those coordinates. You
can get the standard errors of the predicted value as well:

> predict(topo.lo, data.frame(x = 4.8, y = 3.1), se = TRUE)
$fit
[1] 824.0046

$se.fit
[1] 7.677035

$residual.scale
[1] 18.95324

$df
[1] 34.00484

And of course, you aren't restricted to doing this one point at a time:

> predict(topo.lo, data.frame(x = c(4.8, 4.9, 3.1, 2.6),
+                             y = c(3.1, 2.3, 4.5, 5.6)),
+         se = TRUE)
$fit
[1] 824.0046 849.2514 760.2926 744.2987

$se.fit
[1] 7.677035 7.127979 6.364493 7.093619

$residual.scale
[1] 18.95324

$df
[1] 34.00484

Is this what you were looking for?

HTH

G
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Gavin Simpson                 [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%