[R-sig-Geo] GSTAT - can a spline be used in the drift?

[Edzer Pebesma]

Dear Zev, thanks for reminding me about his (off-line).

I looked into this issue, and found that predict.ns (the predict method
that comes with ns usage in predict mode) simply calls ns, and so
re-computes the whole spline with predict data ranges. As the spline
fitted depends on the data that is fed to it, the solution depends on
the data fed to the predict statement, and I believe that's what you noted.

The difference in the solutions seems to be a consequence of the
placement of the knots. The knots are placed along quantiles of the x
values (dist), so if the range differs from one call to the other, the
solution will change.

I'm not sure how bad all this is. The variogram fitted uses a different
x range than the predict step, so the residuals are not identical. Also,
it's a bit muddly, as you're combining two basically different sets of
splines, one based on data, one on prediction values. As long as
quantiles for both sets are identical, results should be identical. Once
their quantiles start being very differently, the differences will increase.

The solution I could find is to fix the knots in calls to ns, e.g. by
the example below.

Hth, with best regards,
--
Edzer

library(splines)
data(meuse)
data(meuse.grid)

coordinates(meuse) = ~x + y
coordinates(meuse.grid) = ~x + y

lznr.vgm = variogram(log(zinc) ~ ns(dist,3), meuse)
lznr.fit = fit.variogram(lznr.vgm, model = vgm(1, "Exp", 300,+ 1))
plot(lznr.vgm, lznr.fit)
Bk = range(meuse.grid$dist)
k = quantile(meuse.grid$dist, 0.5)
preds<-krige(log(zinc) ~ ns(dist,Boundary.knots=Bk,knots=k), meuse,
meuse.grid, lznr.fit)
preds["var1.pred"][1:5,] # these predictions will be different from
those below

data(meuse.grid)
meuse.grid<-meuse.grid[1:5,]
coordinates(meuse.grid) = ~x + y
preds<-krige(log(zinc) ~ ns(dist,Boundary.knots=Bk,knots=k), meuse,
meuse.grid, lznr.fit)
preds["var1.pred"][1:5,] # these predictions are different from those above



On 03/15/2010 09:46 PM, Zev Ross wrote:
> Hi All,
> 
> Sorry for the re-post, I didn't get an answer the first time around and
> thought I'd give it another try. I am noticing in GSTAT that if I
> include a spline as part of the drift in an external drift kriging model
> that the prediction, oddly, depends on how many observations I have in
> the "newdata" argument.
> 
> Using the meuse data as an example, you'll see that the predictions at
> the same first 5 records will be different if I include just those five
> records (second example below) or a larger dataset that includes those
> five records.
> 
> Perhaps using a spline in the drift is not supported and I should use a
> polynomial? Can anyone explain?
> 
> Zev
> 
> 
> library(splines)
> data(meuse)
> data(meuse.grid)
> 
> coordinates(meuse) = ~x + y
> coordinates(meuse.grid) = ~x + y
> 
> lznr.vgm = variogram(log(zinc) ~ ns(dist,3), meuse)
> lznr.fit = fit.variogram(lznr.vgm, model = vgm(1, "Exp", 300,+ 1))
> plot(lznr.vgm, lznr.fit)
> preds<-krige(log(zinc) ~ ns(dist,3), meuse, meuse.grid, lznr.fit)
> preds[["var1.pred"]][1:5] # these predictions will be different from
> those below
> 
> data(meuse.grid)
> meuse.grid<-meuse.grid[1:5,]
> coordinates(meuse.grid) = ~x + y
> preds<-krige(log(zinc) ~ ns(dist,3), meuse, meuse.grid, lznr.fit)
> preds[["var1.pred"]][1:5] # these predictions are different from those
> above
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de