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.pebe...@wwu.de _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
